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INTRODUCTION 


This  report  documents  the  mathematical  details  of  a  solution  to  the 
following  problem:  given  a  general  surface  of  revolution,  produce  precision 
information  off-line  that  a  rotating-traversing  (2-axis)  filament  winding 
machine  needs  in  order  to  wrap  (cover)  the  surface  uniformly. 

Off-line  programming  refers  to  the  use  of  analytical  and  computational 
methods  instead  of  manual  methods  for  teaching  the  machine  the  motions  it  needs 
to  perform  in  order  to  uniformly  cover  a  surface  with  wrapped  filament.  Off¬ 
line  programming  is  important  because  winding  data  can  be  produced  almost  imme¬ 
diately  upon  specification  of  the  surface  profile  function.  Also,  the  filament 
winding  angle  can  be  tailored  to  some  extent  and  various  different  paths  may  be 
tested  on  various  different  surfaces  without  requiring  the  production  of  the 
physical  surface. 

On-line  programming,  on  the  other  hand,  almost  completely  avoids  any  need 
for  mathematical  analysis,  but  it  does  so  at  a  price.  Manual  teaching  is  time 
consuming,  tedious,  error  prone,  less  flexible,  and  it  removes  a  machine  from 
production  while  it  is  being  manually  put  through  its  paces.  Also,  the  very 
first  requirement  of  manual  teaching  is  the  production  of  the  physical  surface 
or  mandrel . 

Various  important  mathematical  areas  involved  in  this  report  are  differen¬ 
tial  geometry,  variational  calculus,  asymptotic  expansions,  numerical  quadra¬ 
ture,  error  analysis,  finite  differences,  numerical  approximation,  number 
theory,  and  of  course,  computer  programming. 

GEODESIC  MINDING  PATHS 

The  curve  on  the  surface  along  which  we  wish  to  wind  a  composite  filament 
will  be  referred  to  as  a  winding  path.  The  ideal  winding  path  is  a  geodesic 
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because  the  filament  will  have  absolutely  no  tendency  to  slip  subsequent  to  its 
meeting  with  the  surface.  If  there  is  no  friction  between  filament  and  surface, 
a  purely  geodesic  winding  path  is  mandatory.  If  there  is.  considerable  friction 
between  filament  and  surface,  one  may  be  able  to  deviate  considerably  from  a 
geodesic  winding  path  with  no  ill  effects.  In  any  case,  the  geodesic  is  the 
logical  starting  point. 

A  good  operational  definition  of  a  geoaesic  path  is  simply  the  path  of 
least  length  between  two  (nearby)  points  on  the  surface.  A  geodesic  in  the 
plane  is  a  straight  line.  A  geodesic  on  a  cylinder  is  a  helix.  A  geodesic  on  a 
sphere  is  a  great  circle. 

A  surface  of  revolution  is  defined  by  its  profile  or  radius  function  r(x), 
where  x  measures  the  distance  along  the  axis  of  the  surface. 

Using  a  right-handed  (x,y,z)  coordinate  system,  the  equation  of  the  axisym- 
metric  surface  of  revolution  is  given  by 

y*  ♦  z*  ■  r(x)« 

For  any  path  on  the  surface  (not  necessarily  geodesic), 

y  ■  r(x)  sin  0 
z  •  r(x)  cos  9 

where  0  is  the  angular  amount  of  wrap  of  the  path  around  the  axis  of  the  sur- 

. .  * 

face,  beginning  at  the  point  (0,0, r(0})  (0  ■  0). 

The  differential  of  arc  length  for  any  point  on  the  path  is  given  by 

(ds)*  *  (dx)4  ♦  (dy)*  ♦  (dz)4 
but 

dy  »  dr  sin  0  +  r  cos  0  d0 
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and 


dz  ■  dr  cos  0  -  r  sin  0  d0 

Therefore 

(ds)*  ®  (Jx)*  +  (dr  sin  0  ♦  r  cos  0  d0)*  +  (dr  cor,  0  -  r  sin  0  d0)* 
Squaring  and  simplifying,  we  have 

(ds)*  «  (dx)*  ♦  (dr)*  ♦  r*(d0)* 
or 

(ds)*  -  (dx)*  +  (~)*(dx)*  ♦  r*(^)*(dx)*  *  ( 1+r ' *+r*0'  * )  (dx) * 


or 


ds  *  (1+r ' *+r*0* * )^  dx 

Therefore,  in  terms  of  the  behavior  of  r  and  0  as  functions  cf  x,  we  may  calcu¬ 
late  the  length  of  any  path  on  the  x  interval  (a,b)  by 


,b  v 

s  «  /  (1+r' *-*r*0' *  )*  dx 

a 

We  will  not  actually  have  occasion  to  calculate  s,  but  this  integral  for  s  is 
important  in  obtaining  the  0  behavior  of  a  geodesic  path  on  the  surface. 

The  fundamental  problem  of  the  calculus  of  variations  is  to  find  0(x)  which 
minimizes  the  integral 

/b  F(x,0,0')dx 

a 


This  minimization  problem  is  equivalent  to  solving  tha  Euler  differential 
equation 

3F  d  3F 
30  "  dx  307 


but  in  our  case, 


F(x,0,0')  «  (l+r'*+r«0'*)31 
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and  this  expression  is  not  dependent  on  8.  Therefore 


3F  Q  d  3F 
38  "  dx  38' 


and 


Now 


3F 

38 


27  *  c  ■  constant 


3F 

ae" 


r*8,(l+r'*+r*8'*)",<  ■  ♦  r*)"5*  ■  c 


Geodesic  curves  which  pass  through  the  axis  exactly  twice  (at  the  poles)  will  be 
called  meridians.  Nongeodesic  curves  having  a  constant  x  coordinate  will  be 
referred  to  as  parallels  or  hoops. 

Assume  the  geodesic  is  tangent  to  a  parallel  at  some  point  xQ  (the  usual 
starting  point  will  be  x0  *  0).  Such  a  point  represents  a  stationary  point  of  x 
as  a  function  of  8.  Therefore 


dx 

d8 


* 


0  at  x  >  x0 


or 


8’  ■  »  at  x  »  x- 


We  may  therefore  conclude  from 


rM's? r  ♦  r*)'*  « 


that 


Therefore 


c  r(x0)  =  rQ 


1 f r 1 S  .V 

*  r»)  *  «  rc 
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The  point  x0  will  be  called  a  turning  point  of  the  geodesic  and  the 
corresponding  radius  r0  will  be  called  the  polar  radius.  The  turning  points  of 
the  geodesic  represent  the  boundaries  between  the  wrapped  region  of  the  surface 
and  the  unwrapped  regions  of  the  surface.  Given  a  surface  of  revolution,  we  are 
theoretically  at  liberty  to  specify  any  point  we  wish  as  a  turning  point  of  a 
geodesic;  having  done  this,  we  have  determined  the  entire  geodesic  and  the  other 
turning  point  as  well.  These  facts  imply  that  in  order  to  wrap  or  wind  alorg  a 
pure  geodesic  path,  we  must  wind  between  two  precisely  equal  polar  radii  and 
r ( x )  may  not  drop  below  the  polar  radius  r0  in  the  wrapped  region.  This 
situation  is  rather  restrictive,  and  we  will  discuss  how  to  ove-come  it  later. 

Solving 


r*( 


Hr'* 

~9~*~ 


for  9 ' .  we  have 


o 

Therefore,  9  for  a  geodesic  with  polar  radius  r0  is  given  by 

0(x)  «  /*  ----  (liCiiSli)5*  dt 
11  Jo  r( t)  V(t)»-r«' 

This  integral  is  improper,  of  course,  because  the  iivegrand  becomes  infinite  at 
the  turning  points.  Intuition  assures  us,  however,  that  the  integral  does 
indeed  exist  and  we  will  subsequently  develop  an  asymptotic  approximation  for  it 
for  small  x.  9(x)  is,  in  fact,  proportional  to  /x  for  small  x. 


PATH-ANGLE  RELATIONS 

Let  the  hoop  angle  h  be  defined  as  the  angle  between  a  hoop  and  the  path. 
The  angle  between  the  path  and  a  meridian  will  be  called  the  winding  angle  w. 
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More  precisely,  th*.  hoop  angle  h  is  the  angle  between  a  vector  tangent  to  the 
path  and  a  vector  tangent  to  the  hoop. 


Letting  i,  j,  and  k  be  the  usual  coordinate  basis  vectors,  the  unit  tangent 
vector  to  any  curve  on  the  surface  is  given  by 

Tp  «  (idx  ♦  jdy  +  kdz)/ds 

■  (idx  ♦  j(dr  sin  0  +  r  cos  0  dO)  ♦  k(dr  cos  8  -  r  sin  8  d0))/ds 
for  the  particular  case  of  a  hoop  or  parallel,  dx  »  0  *  dr,  therefore 

^h  *  Ur  cos  0  dd  -  kr  sin  8  d0)/ds 
For  a  hoop  we  also  have  ds  ■  rd0,  therefore, 

Th  *  j  cos  0  -  k  sin  0 

Taking  the  dot  or  inner  product  of  these  two  unit  tangent  vectors,  we  have 

cos  h  *  Tp  •  Th 

*  {(dr  sin  0  ♦  r  cos  8  d0)cos  0  -  (dr  cos  0  -  r  sin  0  d0)sin  9}/ds 


■  r0' 


dx 

ds 


*  r0' (1+r' *+r*0' 


2 


) 


% 


■  r*  ( 


1+r'* 

■§T*~ 


♦ 


r*)~Vr 


Therefore 


cos  h  = 


r 


Aside  from  its  simplicity,  the  interesting  thing  about  this  relation  is 
that  it  not  only  holds  when  r  is  not  differentiable,  but  it  also  holds  even  when 
r  is  not  continuous.  This  means  that  if  x  >  a  is  a  discontinuity  of  r,  and  r  is 
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set  equal  to  any  value  between  r(a-o)  and  r(a*o),  the  cosine  of  the 
correspr  »  h  will  still  be  given  by  r0/r.  Since  meridians  and  parallels 
inter'  ''gonally,  the  winding  angle  w  is  given  by 

ro 

sin  w  ■  -- 
r 

These  formulas  for  hoop  and  winding  angles  will  hold  (in  slightly  modified  form) 
for  subsequently  defined  nongeodesic  paths  as  well. 

QUASI -GEODESIC  PATHS 

The  following  surface  is  guaranteed  to  be  wrappable  on  a  purely  geodesic 
path.  We  are  considering  just  one  circuit  here,  beginning  at  the  left-hand 
turning  point  (x«0),  winding  to  the  right-hand  turning  point  (x-L),  and 
returning  to  the  left-Nand  turning  point. 


Figure  1.  Equal  polar  radii  for  pure  geodesic. 

The  reasons  why  this  surface  is  guaranteed  .to  be  geodesically  wrappable 
are 

1.  r(x)  ■  r0  at  exactly  two  points  (x*0  and  x=L) 

2.  r(x)  >  r0  for  0  <  x  <  L 

3.  r"(x)  <  0  for  0  <  x  <  L 
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We  will  address  the  problem  presented  by  violating  condition  #1  in  this 
section,  leaving  the  problem  of  violating  condition  #3  for  the  next  section. 

Suppose  we  wished  to  wrap  a  surface  similar  to  the  previous  one,  but  we 
wanted  two  unequal  polar  radii  rQj  and  rg2  (rg2  >  rpi ) * 


wrapped 


Figure  2.  Overwrap  for  unequal  polar  radii. 

If  we  want  rgj  and  rg2  to  correspond  to  turning  points,  the  task  is  mathemati¬ 
cally  impossible  if  we  insist  on  a  purely  geodesic  path.  One  alternative  we 
have  in  this  case,  however,  is  to  wind  a  pure  geodesic  farther  out  on  the 
right-hand  end  of  the  surface  than  we  want  to  and  subsequently  cut  off  the 
unwanted  surface  to  get  our  larger  right-hand  polar  radius.  This  alternative  is 
not  generally  viable,  however,  because  if  the  resulting  composite  surface  is  to 
be  used  as  a  pressure  vessel,  the  cutting  of  the  fibers  or  filaments  may  not  be 
desirable  for  reasons  of  structural  integrity. 

There  is  yet  another  alternative,  however,  which  gives  us  considerable 
flexibility  in  design  and  fabrication  of  the  composite  surface.  We  use  the 
simple  but  elegant  expedient  of  thinking  of  rQ  not  as  merely  a  number,  but  as  a 
function  rQ(x).  Our  quasi -geodesic  path  is  therefore  defined  by 


0'(x) 


r0(x) 

*  r7x7~ 


r (x) *  -  r0(x)*' 
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where  r(0)  »  ro(0),  r(L)  »  r0{ t),  and  r{x)  >  r0(x)  for  0  <  x  <  L.  The  invented 
function  r0(x)  will  be  referred  to  as  the  polar  radius  function. 

The  schematic  for  the  previous  surface  could  therefore  look  like  Figure  3. 


Figure  3.  Quasi -geodesic  with  unequal  polar  radii. 

Here,  we  have  our  turning  points  exactly  where  we  want  them  and  our  polar  diame¬ 
ters  are  exactly  what  we  want  them  to  be  By  using  the  linearly  varying  polar 
radius  function  r0(x),  we  have  given  up  one  small  thing  -  a  pure  geodesic  path. 
Our  path  is  now  geodesic  nowhere,  but  very  nearly  geodesic  everywhere !  By 
inventing  the  polar  radius  function,  we  have  in  one  fell  swoop  generated  an 
infinite  number  of  similar  geodesics  and  selected  a  single  point  from  each! 


We  might  also  have  defined  r0(x)  as  shown  in  Figure  4. 
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Here,  we  have  allowed  r0{x)  to  be  constant  near  the  ends  of  the  surface;  when¬ 
ever  r0(x)  is  constant,  a  pure  geodesic  wrap  results. 

If  we  wanted  greater  hoop  strength  in  the  resulting  composite  surface,  we 
could  define  r0(x)  as  shown  in  Figure  5. 


wrapped 


Figure  5.  Polar  radius  function  for  higher  hoop  strength. 

If  we  wanted  greater  axial  stiffness,  we  might  define  r0(x)  as  shown  in 
Figure  6. 


wrapped 


Figure  6.  Polar  radius  function  for  higher  longitudinal  stiffness. 
Needless  to  say,  the  further  r0  deviates  from  a  constant  function,  the 
further  the  path  will  deviate  from  a  pure  geodesic.  Too  abrupt  a  deviation  from 
a  geodesic  path  can  produce  slipping  and  catching  of  the  filament,  producing  a 
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nonuniform  surface.  For  a  convex  surface,  as  pictured  in  the  previous  drawings, 
slippage  is  the  only  problem  we  can  have  with  a  quasi -geodesic.  For  a  nonconvex 
surface,  another  problem  can  arise  -  lift-off  or  bridging,  where  the  filament, 
in  some  region,  lifts  off  the  surface  and  forms  a  bridge  between  two  points  on 
the  surface. 

The  following  surface  could  suffer  from  bridging  in  the  region  indicated: 


wrapped 


Figure  7.  Nonconvex  surface,  possible  bridging. 

The  avoidance  of  bridging  will  be  the  subject  of  the  next  section. 

WRAPPABLE  PATHS  AND  SURFACES 

As  mentioned  in  the  previous  section,  all  convex  surfaces  are  wrappable, 
provided  we  select  r0  properly  and  no  slippage  problems  occur.  When  r”(x)  >  0 
in  any  region,  however,  we  may  get  bridging  -  where  the  filament  wants  to  lift 
off  the  surface  and  form  a  bridge  between  two  points  on  the  surface.  If 
bridging  can  occur,  the  path,  geodesic  though  it  may  be,  it  said  to  be  unwrap- 
pable.  In  order  to  avoid  bridging,  it  is  first  necessary  to  develop  a  mathemat¬ 
ical  condition  for  the  occurrence  of  bridging. 
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The  gradient  vector  to  the  surface 

S(x,y,z)  ■  y*  +  z*  -  r(x)*  =  0 

will  always  point  perpendicularly  outward  from  the  surface.  The  path  curvature 
vector  may  point  either  inward  or  outward,  however.  As  the  curvature  vector  K 
approaches  the  tangent  plane  from  below,  we  are  approaching  a  bridging 
situation.  If  the  curvature  vector  subsequently  points  outward  from  the  tangent 
plane  at  any  point,  we  have  a  definite  bridging  situation.  The  condition  which 
must  t^refore  hold  in  order  to  avoid  bridging  is  that  the  inner  or  dot  product 
of  the  gradient  vector  and  the  curvature  vector  must  remain  negative. 

A  vector  proportional  to  the  gradient  of  surface  S  is  given  by 

G  ■  -  irr'  ♦  jy  +  kz 

If  P  *  ix  +  jy  ♦  kz  is  a  point  on  the  winding  path,  the  curvature  vector  of  the 


path  is  given  by 


K 


ds*  J  ds* 


♦ 


Therefore 


G*K 


rr.  *i*  ♦  v  lit  +  z 

ds*  y  ds*  z 


d*z 

ds* 
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Using  y  »  r  sin  9  and  z  *  r  cos  0,  doing  the  differentiations  and  simplifying, 
we  find  that 


rr'  -  yy'  -  zz'  *  0 
yy "  ♦  zz"  ■  rr"  -  r*9'* 

Therefore 

G«K  ■  r(r"-r0' * )/s' * 

Since  only  the  sign  of  this  function  matters,  we  define  the  bridging  function  to 

be 


B(x)  ■  r"  -  r0'* 

where  bridging  takes  place  if  B(x)  >  0  and  not  if  B(x)  <  0. 

While  it  is  already  intuitively  obvious  that  a  convex  surface  cannot 
experience  bridging,  it  is  now  also  mathematically  obvious  because  B  cannot  be 
positive  if  r"  is  negative.  On  the  other  hand,  if  the  surface  is  concave  (r"  > 
0)  in  any  region,  B  can  be  positive  if  0'  is  not  sufficiently  large.  In  order 
to  avoid  bridging,  we  simply  make  sure  that  0'  will  always  be  sufficiently  large 
by  making  rQ  sufficiently  large. 

Writing  the  condition  for  no  bridging: 

r"  -  r0'*  <  0 

and  inserting  the  expression  for  0',  we  have 


r" 


r  • 


-2  .  —licll  <  0 

r*  r*  -  r*  '  u 
0 


Isolating  r0  on  one  side,  we  have 


r*r"  v 
ro  >  r>> 
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We  define  the  right-hand  side  of  this  inequality  to  be  the  polai  radius 


lower  bound  function  A(x) 


A<x) 


0  if  r"  <  0 


( _ 

'l+r' *  ♦  rr,,r 


if  r"  >  0 


We  see  therefore,  that  bridging  can  be  avoided  if  we  make  certain  that  r 
A(x )  everywhere.  Note  that  if  r"  ■  ®  (a  sharp  valley),  then  A  *  r. 

PATH  BEHAVIOR  NEAR  TURNING  POINTS 


In  this  section  we  obtain  one-term  asymptotic  expressions  for  9'(x) 
vicinity  of  turning  points  (<?'  «  »).  The  general  expression  for  9'(x)  i 
by 


e'(x) 


.2!?!  /..lidixli _ 

r(x)  V*(x)  -  r* (x) ' 


For  x  near  zero  we  have 

r(x)  ~  r(o)  +  xr’ (o) 
and 

r0(x)  -  rc{o)  ♦  xr^  (o) 

Now 


r(x)»  -  r0(x) *  =  (r(x)-r0(x))(r(x)+r0(x)) 
but 


Therefore,  we  have 


r0(o)  »  r(o) 


r(x)  -  r0(x)  -  x(r‘ (o)-r^(o) ) 

and 

r(x)»  -  r0(x) *  -  x(r' (o)-r^(o  '(2r(o)) 


o  ( x )  ^ 


in  the 
s  given 
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Ultimately, 


9*{x) 


•_*o(o)  . _ _ 

*r(o)  2xr(oJ{r”io)-r^IoJ ) 


and 

For  x  near  L  we  have 

r(x)  -  r(L)  ♦  (x-L)r'(L) 
and 

r0(x)  ~  r0(L)  ♦  (x-L)r^(L) 
but 

r0(L)  *  r(L) 

Therefore 

r(x)  -  r0(x)  ~  (x-L) (r ' (L)-r^(L) ) 
and 

r(x)*  -  r0(x)*  -  (x-L) (r' (L)-r^(L) ) (2r(L) ) 

-  2(L-x)r(L)(r;(L)-r'(L)) 

and  ultimately. 


.1  ( . Li_c  :i?i! 

ft  V2r (o) (.  1 (o)-r^(o) ) ; 


e(x)  - 

r(o) (r 1 (o)-r' (o) ) ' 


fl'(x) 


ro(L)  .  1  ♦  r'XL]* 

r(L)  '2(L-x)r7L)(r^(L)-r”(L)) 


,  ( - _ 

'2r(L)  (rg(L)-r'  (L) ) ' 
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We  see,  therefore,  that  9'(x)  has  square  root  singularities  at  the  left  and 
right  turning  points.  Integration  of  square  root  singular  functions  will  be 
discussed  in  the  next  section. 

A  SQUARE  ROOT  SINGULARITY  QUADRATURE  FORMULA 

In  this  section  we  will  develop  a  one-point  Gaussian  quadrature  formula 
with  error  term  for  integrating  square  root  singularity  functions  in  the  vicin¬ 
ity  of  the  singularity.  We  have  the  following  integral  to  evaluate: 

I  ,  /a+h  fill  dt  for  a  >  0 
a  /~t 

We  rewrite  this  integral  as  two  separate  integrals: 

I  .  /*  fill  dt  .  ,*  fill  dt 

a  /t  a+h  /t 

We  will  use  the  following  formula  for  integration-by-parts  on  numerous  occa¬ 
sions: 

/b  f (x)g(x)ux  *  f(b)/&  g(x)dx  -  /b  f'(x)/*  g(t)dtdx 
a  a  a  a 

Using  integration-by-parts  on  I,  we  have 

I  *  f (x )  /  f^dt  -  /  f'(t )/  u~^dudt 

a  a  a 

-  {f (x)/  t-5*dt  -  /  f'(t )/  u'^duat} 
a+h  a+h  a+h 

■  2f  (x )  (x’i-a5*)  -  /X  2f '  ( t )  (t^-a^)dt 
a 

-2f(x)(x5Ha+h)*)  ♦  /X  2f'(t)(t*-(a+h),i)dt 
■  a+h 
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2f (x)  ( (a+hl’f-a5*) 


Let 


Therefore 


-  /*  2f  {t)(tJ*-a5i)dt 

a 

+  /'  2f 1 (t) (t^-{a+h)^)dt 
a+h 


E  «  I  -  2f(x)((ah)%-a*) 


-  I  - 


__2hf]x}___ 
(a+h)5*  +  a* 


f  *  ■■  J*  f '  (tHt^-a^Jdt 

**  of 


♦  /X  f'(t)(t*-(a+h)*)dt 
a+h 

Using  integration-by-parts  again  on  the  last  two  integrals  gives  us 

|  ■  -  |f'(x)/X  t5*  -  a^dt  -  / *  f”(t)/t  u^  -  a^dudt) 

L  a  a  a 

♦  f * (x ) /  tH  -  (a+hj^dt  -  /  f"(t)/*  -  (a+h^dudt 

a+h  a+h  e+h 

«  -  f,(x)(£/jt*/*-aJJt)|  X  +  /X  .M(t)(*/»u»/*-^u)|  dt 

a  a  a 

XX  t 

+  f’(x)(a/9t»/*-(a+i.)*t)|  -  /  f"{t)(*/?u»/*-(a+h)^u){  dt 

a+h  a+h  a+h 

«  -f  ’  (x)  (*/sx3/*-aJix-a/sa3/ *+a3/*) 


x 

+  /  f"(t)  (*/»t3/*-aJit-*/3a3/*Ta3/*)dt 

a 

+  f,(x)(*/3x3/*-{a+h)%x-*/3(a+h)3/*  +  (a+h)3/*) 

+  /a*h  f"(t)(,/3t3/*-(a+i  ,'Jft-j/s(a+h)3/*  +  (a+h)3/*)dt 
x 


l*i 


Therefore 


|  *  flxllxla^-la+h)^  +  1/3 ( (a+h)3/*-a3/* ) ) 

♦  /X  f"(t)(i/st3/*-a*t+i/*a3/*)dt 
a 

,a+h  ,  - 

+  /  f"{t)  (2/3W  *-(a+h)^t+i/s(a+h)9/*)dt 

x 

If  we  now  pick  x  to  zero  the  f'(x)  term,  we  will  get  a  quadrature  formula  which 
is  exact  for  linear  f.  Therefore,  let 

xfa’Ha+h}*)  +  i/s( (a+h)s/*-as/ * )  *  0 

Solving  for  x  gives 

3<(a+h)*-a*)  "  "  ~~3(a*h-a) 

{a+h)*  -  a3/*{a+hl*  ♦  a^a+hl^/^a* 

"  .  ~3h . 

2ah+h*  +  a^Ja+h^Ja+h-aJi 
. 3h . 

■  (2a+h+(a*+ah)*)/3 

We  therefore  have 


where 


and 


, .  rh  mi  dt ,  .ffiUsi  ,  e 

a  (a+h)<+a< 

x  *  i/a(2a+h+(a*+ah)Jf) 

I  *  /*  f"(t)  ( */3ts/*-a5it+i/aa3/*  )dt 
c  a 

♦  /9+h  f"(t){i/3t3/*-(a+h),it+i/3(a+h)3/*)dt 
x 
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We  note  here  that  if  a  =  0,  the"  x  *  ^  ,  and  for  large  a  (relative  to  h), 

(a+h)^  *  a^d  ♦  ~  a’id  ♦  ~) 

H  jnce 

x  -  i/a (2a+h+a(l  ♦  |~)) 

h 

*  i/*(2a+h+a+  ^ / 

.  3h.  h 

*  t/s(3a  ♦  2~)  *  a  ♦  2 

We  see,  therefore,  that  for  large  a,  this  quadrature  formula  becomes  the  mid¬ 
point  rule. 

We  now  devote  the  remainder  of  this  section  to  developing  an  upper  bound  on 
the  absolute  error  |  E  |  .  From  the  expression  for  E,  define: 

♦  (t)  *  i/st*/*  -  a^t  ♦  i/sa»/ * 
and 

*( t)  »  */3tJ/*  -  (a+h)%t  ♦  i/s (a+h) »/* 

We  first  prove  that  <p  and  v  are  non-negative  in  their  respective  ranees  of 
integration.  This  is  important  when  we  apply  the  absolute  value  triangle  ine¬ 
qualities  to  the  expression  for  E. 

♦(a)  *  t/3 a*/*  -  a*/*  ♦  *  0 

<fc'  (t)  ■  tH  -  a5*  ,  ip'  (a)  *  0 

$"(t)  ■---><) 

2/t 

Hence,  <p  and  <p'  are  zero  at  a  and  $  is  concave  upward  to  the  right  of  a.  Tliis 
shows  that  ♦  is  positive  to  the  right  of  a. 
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<Ma+h)  =  */i(a+h)3/ 1  -  (a+h)3/*  ♦  i/s(a+h)3/*  *  0 

^'(t)  =  t5*  -  {a+hjH  ,  i(»'(a+h)  =  0 

U>"(t)  *  >  0 

2/t 

Hence,  t|>  is  positive  to  the  left  of  a+h. 

Using  the  absolute  value  triangle  inequalities  for  sums  and  integrals  on 
the  expression  for  |  gives  us: 

I  |  M  /X  f"(t)<Mt)dt  ♦  /a+h  f"(t)<*(t)dt  | 

<  I  r  f"(t)*(t)dt  !  +  1  /a+h  f-{t)i|r(t)dt  I 

a  x 

<  /X  I  fn(t)  |  ♦(t)dt  +  /8+h  |  f»(t)  |  4>(t)dt 

a  x 

<  If  "I  (/*  <Mt)dt  +  /8+h  tp(t)dt) 
a  x 


where 


If ’* *  *  sup{|  f"(t)|  :  a  <  t  <  a+h} 


Calculating  the  integrals  of  $  and  4»: 


X  X 

/  4»(t)dt  =  /  */3ts/*  -  a*t  +  i/sa*/*  dt 

a  a 

*  4/i*(x*/*-a*/*)  -  \a%(x*-a*)  +  i/sa3/*(x-a) 
a+h  a+h 

/  <l»(t)dt  =  /  */st3/*  -  ( a+h J^t  +  i/s(a+h)3/ *  dt 

x  x 

4/n((a+h)*/*-x*/*)  -  ^(a+h)H( (a+h)*-x*)  +  i/3 (a+h) */* (a+h-x) 


We  therefore  have 


|||<  lf"l{ «/i * ( (a+h) ,/*-a‘/* )  -  JjxMa’Ha+h)*) 
♦  i/sx(a3/*-(a+h)3/*)  ♦  \a*/*  -  ^/»a*/^ 

-  %(a+h)»/*  +  i/s(a+h)*/*{ 


Therefore 


15  j  E  |  <  lf"B{8((a+h)*/*-a»/*) 

-15x*{a*-(a+h)*)  ♦  i0x(aJ/*-(a+h) J/* ) 

+15a*/»  -  10a*/*  -  15(a+h) */*  ♦  10(a+h)*/*J 
»  If  "I  {l5x*  ( (a+hj^-a3*) 

-  10x( (a+h)*/*-aJ/ * )  +  3{a+h) */*-3a*/* } 

Recall  now  that 

x  *  i/sUa+h+afya+h)5*) 

How  let  k  *  a/h,  p  «  k3*,  and  q  =  (k+1)3*.  Therefore,  a5*  *  ph3*  and  (a+h)3*  *  qh3*. 
Hence 

x  *  i/j (2kh+h+ph^qh^) 

■  §  (2k+l+pq) 
and 

x*  a  (4k*+l+p*q*+4k+4kpq+2pq) 

»  (4k*+l+k(k+l)+4k+pq{4k+2) } 

«  (5k*+5k+l+pq(4k+2) ) 

Our  error  bound  therefore  becomes 

15  |  E  |  <  If "I { 15  •  jj-  (5k*+5k+l+pq(4k+2) ) (qh^-ph3*) 

-  10  •  £  (2k+l+pq) (q3hs/*-p3h*/* ) 

♦  3q*h*/*  -  3p*h*/ * } 

Therefore 

45  |  E  |  <  h*/*  If "I { 5 ( q-p ) (5k*+5k+ l+pq(4k+2) ) 

-  10 ( (k+1 )q-kp) (2k+l+pq)  +  9(k+l)*q  -  9k*pj 
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Simplifying  further,  we  have 


Therefore 


where 


and 


But 


Now 


and 


45  |  E  |  <  h*/*  Ilf  "H  {q(5(5k*+5k+l  )-5k(4k+2 ) 

-  10(k+l)(2k+l)  +  10k*+9(k+l ) * ) 

+  p(-5(5k*+5k+l )  +  5(k+l)(4k+2) 

♦  10k(2k+l)  -  10 (k+1 } 2  -  9k* ) J 
■  h*/* If" I {q(25k*+25k+5-20k*-10k 
-  20k*-30k-10+10k*+9k*+18k+9) 

+  p(-25k*-25k-5+20k*+30k 
♦10+20k*+10k-10k*-20K-10-9k* ) } 

45  |  E  |  <  h*/*l|f,,||{Aq  -  Bp} 

! 

i 

A  «  4k*  ♦  3k  +  4 

'  i 

j 

B  *  4k*  +  5k  +  5 

! 

Aq  -  Bp  »  i*9:§Pll*9!!P) 

K  Aq  +  Bp 

i 

Afgf _-_Bf pf  | 

Aq"~Bp 

«  _z_kBf 

”Aq"+~Bp 

A*  =  (4k*+3k+4) *  *  16k4  +  24k1  +  41k*  +  24k  +  16 

B*  =  (4k*+5k+5) *  =  16k4  ♦  40k1  65k*  +  50k  +  25 
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Therefore 


and 


(k+1 )A*  *  16k*  +  40k4  +  65k3  +  65k*  +  40k  +  16 


kB*  ■  16k*  +  40k4 


We  therefore  have 


(k+l)A*  -  kB*  * 
We  can  now  summarize  the  results  of  thi 


+  65k3  ♦  50k*  +  25k 

15k*  ♦  15k  +  16 
s  section: 


where 


and 


where 


till  dt  ,  ,  E 

a  Ji  P  *  q 


(a  )  0) 


k  «  a/h,  p  x  /k,  q  a  /kTI 
x  *  ^  (2k+l+pq) 


4  1/41  h»/*lf'‘Mj--5-) 

Aq+Bp' 


and 


A  *  4k*  ♦  3k  ♦  4 
B  ■  4k*  ♦  5k  ♦  5 
C  x  15k*  ♦  15k  +  16 


An  approximate  bound  on  the  local  relative  error  is  given  by 


Letting 

we  see  that 


<  ..e.i, 9__ 

-  2b5*  |  f  (x 


i/«*h*/*l!f"ll{r----) 

Aq+Bp' 


■  i/»«h* 


If”! 

Hoof 


'Aq+Bp  ' 


Q(k) 


Sle±9l. 

Aq  +  Bp 


Q(o)  x  4  and  Q(»)  *  3.75. 
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\ 


Hence,  R  is  relatively  independent  of  k  and  we  have  approximately 


If"  l 

R  i  »/•**•  pfIJyr 

INTEGRAL  INVERSION 

In  this  section,  we  will  compute  the  inverse  of  the  integral  of  a  positive, 
continuous,  piecewise  linear  function.  In  subsequent  sections,  we  will  have  two 
distinct  applications  for  this  process.  We  begin  with  the  points 

(Xi.fi)  1  <  i  <  n 

where  f^  >  0.  Define  the  continuous,  piecewise  linear  function  f  on  the 
subinterval  such  that 

X  *  X  4 

f  i  <x)  =  fi  +  (fi  +  1-fi)  x-i  <  X  <  X-j+J 

where  h^  *  xi+j  -  x,.  Therefore 

^ i (xi )  =  ^i 

and 

fi<xi+l)  *  fi+l 

The  integral  of  this  function  is  giver,  by 

F(x)  *  /*  f(t)dt 
X1 

Our  goal  is  to  compute  the  inverse  of  F  (F-1),  i.e.,  if  F ( a )  =  b,  then  a  = 
F_1(b).  Now,  if 

xi  <  x  <  Xi+1 
xi  x 

F(x)  *  /  f (t )dt  +  /  f ( t )dt 

Xj  Xi 

*  HXi)  ♦  /*  fi  ♦  (fi+1-fi)dt 

X^  ’*1 


24 


(x-x-j)* 

*  Fi  ♦  (x-Xi)fi  ♦ 

Therefore,  if  we  are  given  F  such  that  F,-  <  F  <  Fi+i,  we  need  only  solve 
the  following  quadratic  equation  for  x  in  order  to  invert  the  integral: 

{fi+l-fi ) 

'"2hi  '  (x'xi>*  +  fi(x-xi)  -  (F-Fi )  >  0 

Therefore 

-fi  ±  (f *^2/hi (f i+!-f , ) (F-Fi ) 
x  -  X-i  ■ - 

*i+l  -  fi 
"hi 

If  f.j+i  >  fi,  we  need  the  plus  sign  in  order  to  get  the  positive  root.  If  f^i 
<  f.j,  we  still  need  the  plus  sign  in  order  vo  get  the  smallest  root.  Therefore 

{✓5-f i )hi 

1  f i+1  -  fi 

where 

0  -  fi  ♦  fiT  (f i+l"f i ) (F-Fi ) 

Noting  that 

Fi+1  -  Fi  •  (fi*fi+i) 


D  can  be  written 


0  -  fi  *  l~.  (fi+l-fi )  (F-Fi ) 


hi/2 ( f i+f i+i )  , 

”Fi:r=“ *  1 


(f*+i-fix 


F-Fi 

F^rF 


T) 

i 


Letting 


clearly 


P 


f i+l  "  Fi 


0  <  p  <  1 
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AUTOMATIC  PIECEWISE  LINEAR  APPROXIMATION 

In  the  process  of  generating  the  data  which  the  filament  winder 
understands,  occasions  arise  when  we  must  take  a  known  (computable)  function  and 
generate  points  on  the  function  which  yield  a  good  piecewise  linear  approxima¬ 
tion  to  the  function.  The  simplest  way  of  doing  this  is  to  generate  a  set  of 
equally-spaced  values  for  the  independent  variable  and  evaluate  the  function  on 
this  uniform  mesh.  This  is  a  poor  method  to  use  here,  however,  because  the 
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functions  we  are  dealing  with  can  have  regions  of  smooth,  regular  behavior 
followed  by  regions  of  abrupt,  irregular  behavior.  A  uniform  mesh  will  evaluate 
the  function  too  often  in  regions  of  low  curvature  and  not  often  enough  in 
regions  cf  high  curvature.  The  resulting  piecewise  linear  approximations  will 
therefore  be  excessively  accurate  vi  regions  of  low  curvature  and  not  accurate 
enough  in  regions  of  high  curvature.  The  problem  of  generating  good  piecewise 
linear  approximations  becomes  even  more  significant  when  there  are  microproc¬ 
essor  hardware  or  software  limitations  on  how  many  data  points  you  may  use  in 
your  digitally  controlled  filament  winder. 

The  technique  used  here  to  generate  good  piecewise  linear  functional 
approximations  is  the  method  of  approximate  error  equidistribution  discussed  in 
general  in  Reference  1. 

The  idea  is  to  create  a  mesh  for  the  independent  variable  (x  here)  which 
approximately  equidistributes  the  error  in  piecewise  linear  interpolation.  V  ■. 
error  bound  for  linear  interpolation  on  the  it^1  subinterval  is 

h* 

e1  """(I) 

where 

hi  ■  *i+1  ~  *i 
and 

lf"l(i)  ■  max  |  f”(x)  | 

X^X^Xi  +  j 

We  would  like  to  select  the  mesh  such  that 

h*lf"l^j  a  constant 


New  York,  1978. 
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or  such  that 


This  is  equivalent  to 


h.j »f "I / .j *  constant 


xi+l 

/  *  constant 


Now,  asymptotically  as  h,  -  0, 

I  ""■(i)*  *  I  f”(x)  i*|-o 

for  any  x  in  the  ith  subinterval. 

We  therefore  determine  the  mesh  by  selecting  the  x's  for  which 

xi+l 

I  |  f”(x)  I  ^dx  a  constant 

X-i 


I  f"(x)  |  *  «  g(x) 


xi*l 

/  g(x)dx  »  C 
xi 


Therefore,  if 


G(x)  a  /  g(t)dt 
X1 


We  therefore  have 


G(x-j )  «  (i-l)c 

G(xn)  *  (n-l)c 

G<xi)  i  -  1 
G(xn)  *  n'-'l 


xi  *  G'M™  G<xn)) 


where  G~’  is  the  inverse  of  G. 
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A  simple  example  is  perhaps  in  order  here.  Suppose  we  wished  to  approxi¬ 
mate  xP  (p  >  0)  by  a  piecewise  linear  function  over  the  interval  (0,xn). 


Therefore 

and  hence 

It  is  easy  to  see  that  if  p  <  2,  the  resulting  mesh  will  tend  to  cluster  more 
points  near  the  left-hand  side  of  the  interval,  while  if  p  >  2,  more  points  will 
cluster  on  the  right.  If  p  *  2,  the  mesh  is  uniform  because  f"(x)  is  constant. 

Usually,  we  do  not  know  f"(x),  because  f  may  only  be  computable  and  not 
given  as  a  simple  explicit  function  of  x.  We  can,  however,  easily  estimate 
f"(x)  over  a  uniform  mesh  using  finite  differences.  We  may,  therefore,  estimate 
g  on  a  uniform  mesh  and  define  our  g  (estimate)  to  be  piecewise  linear.  G  will 
then  be  piecewise  quadratic  and  invertible  by  the  technique  discussed  in  the 
previous  section  on  integral  inversion. 

If  we  let  x,  be  the  i**1  old,  uniform  mesh  point  and  xj<*  be  the  kth  new, 
error  equidistributing  mesh  point,  we  may  use  the  results  on  integral  inversion 
to  define  the  new  mesh: 


f(x)  ■  xP 


f’(x)  oc  x1 
f"(x)  «  x1 


P"1 

P-2 


I  f"(x)  |  *  *  g(x)  *  xp/2_1 


G(x)  oc  /  tP/,2"*dt  oe  Xp/2 


G(xi)  i  .  !  x<  P/2 


G(xn)  n  -  1 


r 

(--) 


i-1  2/P 

xi  “  xn 
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*1*  -  *1 
xn*  a  xn 


where 


xk* 


2(Gk*-Gi) 

xi  ♦ . y 

9i  ♦  VO 


7  <  k  <  n-1 


Gk*  > 


Gn 


P  » 


Gi+1  “  Gi 


r  «  (l-p)g*  ♦  pg*+i 


and  i  is  defined  by 


The  g's  are  defined  by 


where  f^"  is  defined  by 


G-j  <  Gk*  <  G,+i 


9i 


U" 


X 


A  ^ i-1  ”  2f -f  ♦  f -j 

Ti  . h* . 

and  h  is  the  uniform  mesh  spacing. 

The  only  minor  problem  that  can  arise  in  this  process  is  when  some  con¬ 
secutive  values  of  g  are  zero  (and  some  consecutive  values  of  G  are  equal)  in  a 
perfectly  linear  section  of  f.  When  there^ is  at  least  one  nonzero  g,,  however, 
the  problem  of  G^  »  Gk*  ■  G,+i  occurs  with  probability  zero.  We  can  therefore 
avoid  tnis  problem  by  simply  setting  xk*  *  x-j  if  Gk*  *  G-j  and  not  compute  p. 

If  we  have  some  idea  of  how  accurate  we  would  like  our  piecewise  linear 
approximation  over  the  new  mesh  to  be,  we  may  estimate  the  needed  number  of 
points  in  the  new  mesh  in  the  following  manner.  Suppose  we  desire  an  absolute 
error  tolerance  E.  We  would  like 


31 


where 


E  *  |  h*g*  *  constant 
8 


g(x)  -  |  f"(x)  |  * 

as  before. 

In  particular,  we  would  like 

E  *  §  hmin9*ax 

"ft 

Now  there  will  be  an  i  in  the  new  mesh  for  which  xi+1  -  x^  =  hm.jn  will  hold. 
But  we  already  know  that 

G(x*+1)  -  G(x*)  «  G(xn)/(n-l) 

and  that  * 

*  *  ,xi+l 

®(xi+i )  "  «(Xi)  a  /  *  9(x)dx  -  hmingmax 
xi 


Hence,  we  have  approximately  that 

G(xn) 

"min  *  (R:i)g;« 

3i  d 

1  G(xn)«  t 

E  *  §  (n-i)*g*  9ma 

'  '  "max 


Solving  for  n  gives  us 


n  a  1  ♦ 


G(xn) 

TiF 


PATH  COMPUTATION 

In  this  section,  we  compute  the  actual  geodesic  (or  quasi -geodesic)  path  on 
the  surface  of  revolution.  That  is,  we  obtain  the  relation  between  x  and  0. 
Given  the  radius  or  profile  function  r(x)  and  the  polar  radius  function  r0(x), 
we  obtain  9  by  integrating: 
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an  dt 

Ji 


where  a  =  l-x-h  and  f ( t )  *  /t9'(L-t), 


Up  to  this  point,  we  have  been  using  x  as  the  independent  variable  in  this 
analysis.  This  is  the  obvious  choice,  since  the  radius  function  r  and  the  polar 
radius  function  r0  must  be  described  as  functions  of  x.  At  this  point,  however, 
we  will  switch  to  using  9  as  our  independent  variable  because  x  as  a  function  of 
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9  has  no  singularities,  can  therefore  be  subjected  to  higher  order  inter¬ 
polation,  and  is  a  single-valued  function  for  an  entire  circuit  (left  turning 
point  to  right  and  back  to  left  again).  In  order  to  define  x  as  a  function  of  9 
accurately,  we  proceed  with  the  following  steps: 

1.  Generate  a  first  x  mesh  which  is  uniform. 

2.  Integrate  9'  over  this  x  mesh  to  give  the  first  6  mesh. 

3.  Use  this  9,x  data  to  generate  a  second  9  mesh  which  equidistributes 
the  error  in  piecewise  linear  interpolation  of  x. 

4.  Use  piecewise  cubic  interpolation  of  the  first  set  of  9,x  data  and  the 
second  9  mesh  to  generate  a  second  x  mesh. 

5.  Integr  te  9*  again  with  respect  to  the  second  x  mesh  and  produce  the 
third  and  final  9  mesh. 

6.  Define  x  as  a  function  of  9  using  piecewise  cubic  interpolation  of  the 
second  x  mesh  with  respect  to  the  third  9  mesh. 

In  what  follows,  let  dot  (•)  denote  d/d9  and  let  n  be  the  number  of  points 
in  each  mesh.  Ordinarily,  higher  order  interpolation  would  require  a  higher 
order  derivative  than  the  second  for  the  process  of  error  equidistribution,  but 
if  we  define  the  nodal  derivatives  for  the  interpolating  piecewise  cubics  in  the 
following  manner: 

x(e1=o)  =  x(9n)  =  o 


i(0i) 


*i+l  ~_*i-l 
®i+l  "  ®i-l 


(1  <  i  <  n) 


then  the  cubics  will  be  only  0(h2)  accurate  (same  as  linear)  for  a  nonuniform 
mesh,  and  hence,  a  good  mesh  for  piecewise  linear  approximation  will  also  be  a 
good  mesh  for  this  particular  piecewise  cubic  approximation.  Another  reason  for 
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defining  the  nodal  derivatives  this  way  is  to  produce  a  more  stable  (in  the 
sense  of  preserving  monotonicity)  cubic  interpolant  over  nonuniform  meshes. 

The  use  of  a  piecewise  cubic  to  define  x(0)  turns  out  to  be  numerically  critical 
around  the  turning  points;  if  we  do  not  have  x(o)  =  x(0n)  =  0,  we  get  a 
noticeable  kink  in  our  winding  data.  The  piecewise  cubic  interpolant  is  given 
by 


x(0)  *  xi(l-3p*+2p>)  ♦  xi+1(3p*-2p3) 

♦  (fli+l-di)(xi(p-2p*+p3)  +  xi+1(-p*+ps)) 


where 


and 


®i  *  ®  *  ®i+l 

0  -  0, 

P  ®i+l  "  ®i 
*1  »  *n  *  0 


*i*l_:_*i;l 
®i*l  ‘  ®i-l 


(1  <  i  <  n) 


PATH  WINDER  RELATIONS 

The  (2-axis)  filament  winding  machine  or  winder  consists  basically  of  the 
surface,  which  is  rotated  on  its  axis  by  the  winder  and  a  carriage  which  pays 
out  a  taut  filament  from  delivery  point  0,  which  moves  parallel  to  and  at  a  suf¬ 
ficient  distance  from  the  axis  of  the  surface.  The  task  at  hand  is  to  coor¬ 
dinate  the  rotational  movement  of  the  surface  and  translational  carriage 
movement  in  such  a  way  as  to  lay  down  the  filament  on  a  predetermined  path  on 
the  surface.  This  is  accomplished  ultimately  by  computing  carriage  position  as 
a  function  of  surface  rotation.  This  is  the  only  data  that  the  winder 
"understands."  For  purposes  of  visualization,  it  is  beet  to  think  of  the 
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surface  and  the  desired  path  on  it  as  being  stationary  and  the  carriage  as  both 
translating  parallel  to  the  surface  axis  and  rotating  about  the  surface. 

Given  a  point  P  on  the  surface  path,  the  tangent  plane  is  formed  by  the 
vector  tangent  to  the  path  and  the  vector  tangent  to  the  meridian. 


Figure  8.  Tangent  plane  relations. 

The  taut  filament  emanating  from  the  delivery  point  D  just  meets  the  path 
at  point  P.  The  winding  angle  w  is  given  by 


half  circuit)  the  path  point  P.  In  this  picture,  r'  is  assumed  to  be  positive. 
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Figure  9.  Meridian  plane  relations. 


We  immediately  have  the  relations 

c 


*  r*  and  a*  =  I*  ♦  cs 


Therefore 


and 


a*  «  1*  ♦  r ' *1* 
«  I* (1+r » *) 


a*r*  r*i* (l+r ' * ) 

t>*  -  -rr"h' a 

o  o 


The  hoop  plane  is  formed  by  the  hoop  or  parallel  through  P,  The  lead  angle 
A  is  given  by 


cos  A  * 


r  *  c 


where  f  is  the  filament  delivery  height  or  distance  from  the  surface  axis  to  the 
delivery  point  D. 
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Note  that  if  r'  *  0,  1  and  0'  have  the  same  sign.  Therefore, 


1  > 


-rr’  ♦  ✓(fr * ) *  ♦  (f  *-r* )  (r9' )! 

0* 


(r0 '  )*  ♦  r'1 


-rr'  ♦  0'/r*j[f *-r*|  +  {fr'/e']* 

* . (r07 )  *~>~r '  * . 

The  lead  distance  is  therefore  given  by 

*  r*+(r'/e')* 

where 

0  «  r* (f *-r*)  +  (fr'/e*)* 

This  formula  for  I  is  most  accurate  when  r'/0'  <  0.  If  r'/6'  >  0,  we  rational¬ 
ize: 


but 


{DWr'/S'l/e'  05*  ♦  rr'^0' 

*  mT(?V0')5"  *  o*”’rrV0' 

,  ..Ao-Arrwr.nur. _ 

(r*+(r'/0' )* ) (D^+rr '/0) 


D  -  (rr '/0' ) *  «  r*(f*-r*)  ♦  (fr'/0')*  -  (rr'/e1)1 
«  (f *-r* ) (r*+(r'/0’  )* ) 


Therefore 


I  ,  _lf (r*/0‘  >  0) 
rr’/0'  ♦  D* 

We  have  used  the  reciprocal  of  0'  in  these  formulas  because  it  is  always  finite. 
Having  established  I  in  terms  of  the  filament  delivery  height  and  the  path,  we 

i 

may  compute  X  from 


cos  X  ■ 


r  ♦  fr' 
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Finally,  the  carriage  rotation  R  (relative  to  the  delivery  point)  is  given  by 

R  *  0  ♦  A 

and  the  traversing  carriage  position  is  given  by 

T  ■  x  ♦  1 

Subsequently,  it  will  be  necessary  to  express  T  as  an  explicit  function  of  R 
the  only  data  the  winder  "understands"  directly. 


WINOER  OATA  GENERATION 

The  process  of  generating  the  data  that  the  winder  understands  (carriage 
traverse  position  T  as  a  function  of  carriage  rotation  R)  is  complicated  by  the 
fact  that  T  is  not  directly  expressible  or  computable  as  a  function  of  R. 
Instead,  R  and  T  are  both  parametric  functions  of  0,  since  given  0,  x  can  be 
computed  and  given  x,  1  and  A  may  be  computed.  The  solution  to  this  problem  is 
to  find  a  9  mesh  which  will  yield  an  R  mesh  which  will  ultimately  equidistribute 
the  error  in  T  as  a  function  of  R.  We  will  let  dot  (»)  denote  differentiation 
with  respect  to  0  in  what  follows.  Since  our  goal  is  to  equidistribute  the 
error  in  T  as  a  function  of  R,  we  first  write  down  the  approximate  condition 
embodying  this  requirement 

(dR)*  |  |  *  constant  »  8E 

but 

dR  *  Rd0 


We  therefore  want 

(Rd0)*  |  |  =  constant 


but 


dT  dT  d0 
dR  =  d0  dR 


•  • 

T/R 
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and 


&;  ■  i~R  «t/R»  -  &  <t/R)  g  =  <TR-TR)/R> 

Hence 

(Rd0) *  I  (TR-TR)/RS  I  =  (dfl)»  !  T-TR/R  |  *  constant 
Therefore,  our  g  function  for  the  0  mesh  is 

g(0)  -  |  T  -  TR/R  |  * 

Using  this  g  function  in  our  usial  error  equidistributing  process  will  yield  a  0 
mesh  which  will  produce  an  R  mesh  which  will  equidistribute  the  error  in  T  as  a 
function  of  R.  The  major  difficulty  first  seems  to  be  computing  the  derivatives 
in  g,  and  although  this  would  indeed  be  a  cumbersome  task  to  do  analytically,  it 
is  a  simple  task  to  do  numerically.  For  instance,  using  central  differences,  we 
can  obtain  the  following  after  some  simplification: 

T,  -  T^/Ri  *  (2/h«){(Ti  +  1-Ti)(Ri-Pi_1)  -  (T{-T1.1)(Ri  +  1-R1)}/(R1+1-R1.1) 

Me  proceed  algorithmically  as  follows: 

1.  Compute  R  and  T  over  a  dense,  uniform  9  mesh.  For  each  0  value,  com¬ 
pute  x,  1,  and  A  and  store  0+A  and  x+J  in  R  and  T,  respectively. 

2.  Using  finite  differences,  compute  a  piecewise  linear  approximation  to 
9- 

3.  Invert  the  integral  of  g  to  obtain  a  new  0  mesh. 

4.  Evaluate  R  and  T  again  over  this  new  0  mesh. 

To  be  u'able,  the  R/T  data  generated  by  this  process  must  represent  a 
single-valued  function.  This  will  be  the  case  if  the  polar  radius  function  has 
been  properly  defined  to  lie  between  the  radius  function  and  the  polar  radius 
lower  bound  function.  If  the  lower  bound  function  exceeds  the  polar  radius 
function  anywhere  however,  it  is  likely  that  the  Ft  data  will  not  be  monotoni- 
cally  increasing,  causing  unacceptable  loops  or  cusps  in  the  R/T  curve. 
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TIME  BASE  COMPUTATION 


Geometrically  speaking,  the  function  T(R)  is  all  that  is  needed  to  wrap  the 
surface;  but  practically  speaking,  we  must  also  decide  how  R  is  to  behave  as  a 
function  of  time  (or  vice  versa).  Let  us  compute  the  velocity  of  the  delivery 
point  relative  to  the  surface.  The  components  of  the  position  vector  of  the 
delivery  point  are 

x  -  T(R) 
y  *  f  sin  R 
z  ■  f  cos  R 

where  f  is  the  filament  delivery  height. 

The  components  of  velocity  are  therefore 

x  -  T' (R)R 

.  • 
y  a  f  cos  RR 

•  * 
z  «  -f  sin  RR 

where  (•)  denotes  ~  .  The  square  of  the  speed  is 

v*  ■  x*  ♦  y*  ♦  z* 

a  T'(R)*R*  +  f  *R*  cos*  R  ♦  f*R*  sin*  R 
a  ( T * (R) *+f * )R* 

and  the  speed  is 

v  *  (T’(R)*+f*)*R 

If  R  is  a  linear  function  of  t,  R  is  constant  and  |  x  |  will  follow  |  T’(R)  I  •  In 

#  • 
order  to  make  |  x  |  sufficiently  small  everywhere,  we  may  need  to  make  R  so  unac¬ 
ceptably  small  that  we  will  slow  down  the  winding  process  far  too  much  when 
|  T'(R)  |  is  small.  Clearly,  we  need  R  small  when  |  T'(R)  |  is  large  and  vice 
versa.  We  obtain  a  variable  R  by  rewriting  the  previous  equation  as 

dt  a  J  (T'(R)>*f«)4dR 
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Integrating,  we  get 


,^i+l  1  t, 

*1+1  "  ti  *  /R  ■  {T'(R)2+f* )^dR 

Approximating  this  integral  by  the  midpoint  rule  and  using  differences  gives  us 
the  following  time  base  recursion  formula,  where  the  positive  speed  function  v 
is  anything  we  wish: 

*i+l  -  *i  »  ^  ^Ti+l~Ti  >  *  +  f^Ri  +  l-Hi)*)5* 

SURFACE  COVERAGE  RELATIONS 

The  mathematical  machinery  we  have  developed  thus  far  enables  us  to  define 
a  good  quasi -geodesic  path  on  a  surface  of  revolution  and  to  successfully  wind 
filament  along  this  path  by  computing  carriage  traverse  position  as  a  single¬ 
valued  function  of  surface  rotation.  It  is  very  unlikely,  however,  that 
actually  winding  filament  along  such  a  path  will  result  in  a  uniformly  covered 
surface.  Our  task  in  this  section  is  therefore  to  modify  our  nominal  quasi¬ 
geodesic  path  slightly  by  multiplying  0'  by  some  factor  a,  near  unity,  which 
will  guarantee  that  repeated  application  of  the  winding  data  subsequently  pro¬ 
duced  from  the  modified  path  will  serve  to  cover  the  surface  uniformly. 

We  begin  by  noting  the  effect  on  the  polar  radius  function  of  multiplying 
9'  by  o.  The  effective  polar  radius  function  will  be  denoted  by  pQ.  Using  the 
general  definition  of  0,  we  can  write  the  following  equivalence: 

“[o  .1+r^*  m  p_o  .  1+r^*  .5$ 

p  'p2— p*  P  p2-p2 

O  0 

Solving  this  equation  for  pQ  gives  us: 

orr0 

p  - - - - - 

o  (r*+(a*-l)ri)H 
o 
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Note  that  if  a  s  i,  p  s  r0,  and  if  a  =  ®,  pQ  =  p.  It  is  easy  to  prove  that  if  \ 

a  >  1,  then  p0  <  pQ  <  p.  If  a  <  1,  we  may  risk  violating  the  polar  radius  lower 
bound  function.  Remember  that  this  lower  bound  function  is  established  once  and 
for  all  by  the  radius  funct»on  r. 

We  now  define  a  couple  of  discontinuous  functions  for  reference. 

FL(x)  *  Floor  of  x  =  the  greatest  integer  <  x 
CE(x)  =  Ceiling  of  x  *  the  least  integer  >  x  =  -FL(-x) 

Note  that  if  fl(x)  is  the  same  as  FL(x),  but  defined  only  for  positive  x,  then 
the  following  defines  FL(x/  for  all  X: 

FL(x)  =  fl(x+n)  -  n 

whore 

n  a  f  1  (|  x  |  +1) 

We  must  now  investigate  the  closure  properties  of  our  paths.  For  instance, 
if  one  circuit  of  our  nominal  path  wraps  through  a  rational  multiple  of  2tt 
radians,  the  path  is  guaranteed  to  ultimately  close  (i.e.,  return  to  exactly  the 
same  point  on  the  polar  parallel  it  started  at).  If  the  total  wrap  ( 20 ( L ) )  of  a 
nominal  circuit  is  given  by 

"NC  *  2*  O'  ♦  £) 

where  j/k  is  in  lowest  terms,  it  becomes  ODvious  that  the  path  will  repeat 
itself  for  the  first  time  after  exactly  k  circuits.  If  j/k  were  not  in  lowest 
terms,  j  and  k  would  have  a  greatest  common  factor  (GCF)  greater  than  unity. 

Repetition  would  therefore  begin  after  only  k/GCF(j,k)  circuits.  For  the  nomi¬ 
nal  path,  k  could  naturally  be  quite  large.  We  will  establish  uniformly 
wrapping  paths  with  relatively  small  values  of  k.  In  what  follows,  we  will  give 
special  names  and  designations  to  i,  j,  and  k. 
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We  will  call  1  the  number  of  complete  revolutions  per  circuit  and  denote  it 
by  Nrc.  We  will  call  j  the  return  index,  separation  index,  or  number  of  separa¬ 
tions  per  circuit  and  will  denote  it  by  Ng£.  We  will  call  k  the  circuit  index 
or  number  of  circuits  per  repetition  and  will  denote  it  by  NqR.  In  what 
follows,  we  will  also  allow  a  somewhat  looser  interpretation  of  the  term 
"repetition."  For  a  uniformly  wrapping  path,  we  will  say  that  the  first  repeti¬ 
tion  has  occurred  when  the  band  of  filaments  begins  winding  just  alongside  and 
slightly  overlapping  the  first  circuit.  Hence,  a  repetition  may  or  may  not  be 
closed.  We  will  always  use  the  term  "closed"  to  mean  exact  return  to  the  ini¬ 
tial  point. 

The  wrap  of  a  circuit  in  a  closed  repetition  is  given  by 

Ncr 

WCCR  *  <nRC  ♦  fj“") 

The  wrap  for  the  entire  closed  repetition  is,  of  course, 

*CR  *  nCR*wCCR  *  2rr(NCRoNRC+Nsc) 

If  we  actually  wind  a  filament  along  this  path,  the  windings  will  naturally 
close  after  NqR  circuits.  The  pattern  of  filament  bands  on  the  surface  will  be 
such  that  neighboring  left  to  right  travelling  bands  will  be  separated  from  each 
other  by  a  wrap  angle  of  exactly  2ji/NqR.  Neighboring  points  of  tangency  of  the 
bands  to  the  polar  parallels  will  also  be  separated  by  the  same  angle.  We 
therefore  refer  to  2jt/NcR  as  the  separation  angle.  The  angle  2jt  Nsc/Nqr  will  be 
called  the  return  angle  -  the  angle  between  the  band's  point  of  tangency  to  the 
left  polar  hoop  at  the  beginning  of  a  circuit  and  the  band's  point  of  tangency 
to  the  left  polar  hoop  at  the  end  of  a  circuit.  The  return  angle  is  obviously 
the  product  of  the  separation  angle  and  the  separation  index.  We  naturally 
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always  select  values  or  Nqr  snd  Nsc  which  are  relatively  prime,  i.e.,  values  for 
which  GCF  (NgQ,NQp)  *  1.  Another  quantity  which  we  must  introduce  is  the 
number  of  circuits  per  separation,  denoted  by  Nqs.  This  quantity  is  the  number 
of  circuits  which  must  be  wound  in  order  to  obtain  a  net  wrap  advance  of  exactly 
one  separation  angle  beyond  the  starting  point.  Consider  the  following  diagram, 
where  A  is  the  starting  point  on  the  polar  hoop,  B  is  the  return  point  on  the 
polar  hoop  (reached  after  one  circuit),  and  C  is  a  point  on  the  polar  hoop  which 
is  advanced  exactly  one  separation  angle  beyond  point  A. 


Figure  11.  Return  and  separation  angles  for  a  closed  repetition. 

S  is  the  separation  angle  and  R  is  the  return  angle.  Now,  after  NCs^circuits, 
the  path  will  be  tangent  to  the  polar  hoop  at  C.  We  therefore  have 

Ngs,R  *  2nk  ♦  S 
but 


Therefore 


2n_ 

*CR 


NCS*NSC*2,r/NCR  *  2nk  +  2*/Ncr 


4fi 


and 


NCS*NSC  *  1  +  k*NCR 
or 

NCS*NSC  *  1  1,1001  nCR 

Hence,  Ncs  and  Ng^  are  "reciprocals"  in  the  modulo  (or  congruence)  arithmetic 
restricted  to  integers  zerc  through  Ng^-l,  where  multiples  of  NgR  are  equivalent 
(congruent)  to  zero.  We  compute  Ngg  as  the  first  positive  integer  for  which 

(Ncs*NSC-1^NCR 

has  a  zero  remainder.  This  value  of  Ngg  is  guaranteed  to  exist,  provided 
GCF(NSC,NCR)  =  1  (ref  ?). 

We  will  .iow  shew  whai  the  number  Ngg  is  used  for.  Assuming  the  value  of 
N^r  is  relatively  small,  a  single  closed  repetition  will  not  have  covered  much 
of  the  surface  and  we  obviously  won't  cover  any  more  of  the  surface  by  con¬ 
tinuing  to  wind  along  a  closed  path.  We  now  imagine  the  band  of  filaments  to  be 
cut  at  some  point  and  one  cut  end  to  bs  displaced  from  the  other  through  a  small 
wrap  angle  which  we  will  call  the  band  advance  per  repetition  or  Br.  To  get 
uniform  coverage  of  the  surface,  we  would  like  Br  to  be  such  that  an  extension 
of  the  displaced  cut  end  will  continuously  overlap  the  previously  wo^nd  band 
ever  so  slightly.  Suppose  we  were  to  wind  a  couple  more  repetitions,  using  the 
same  Br.  What  would  the  configuration  of  bands  on  the  surface  look  like?  There 
would  still  be  NCR  separated  right  travelling  bands  on  the  surface,  but  they 
would  be  nearly  three  times  as  wide  as  the  first  band.  Also,  there  would  still 
be  two  cut  ends,  separated  by  a  wrap  angle  of  3Br.  Now,  it  should  be  clear  that 
as  the  angular  separation  between  the  two  cut  ends  approaches  the  separation 
angle,  the  surface  will  be  almost  completely  covered  by  a  layer  of  filament.  We 

ZW.  E.  Deskins,  Abstract  Algebra.  MacMillan  Company,  New  York,  1964. 
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will  now  address  the  problem  of  obtaining  exact  closure  of  the  path  after  the 
surface  has  been  covered  by  one  layer  of  filament.  Consider  the  following 
diagram: 


Figure  12.  Starting  points  for  unclosed  repetitions. 

i 

The  points  labeled  one  through  Nrl  (number  of  repetitions  per  layer)  are  the 
starting  points  of  the  first  through  the  Np|_Th  repetitions,  respectively.  Each 

I 

point  is  separated  from  its  neighbors  by  the  band  advance  per  repetition,  BR. 
Now,  for  the  sake  of  uniformity,  we  would  like  to  select  BR  so  that  the  angle 
between  point  C  and  the  point  indexed  by  NRl_  is  exactly  equal  to  Br;  i.e.,  we 
would  like  our  path  to  already  be  winding  directly  on  top  of  itself  when  it 
reaches  point  C  for  the  second  time.  Here,  for  unclosed  repetitions,  the 
separation  angle  S  is  slightly  greater  than  it  is  for  a  closed  repetition.  Here 


where  8c  is  the  band  advance  per  circuit.  We  therefore  want 


♦  NCS*BC  »  Nrl*Br 
but 

Br  «  NCr*Bc 

Therefore 

2tt 

*  Ncs,0c  ■  NRL*NCR*BC 

Hence 

a  . 2  n . 

C  NCR(HCR-NR  -Ncs) 

and 

_  2it 

R  *  NcR"NR[""Nci 

We  have  thus  computed  the  band  advances  necessary  for  our  path  to  reach  point  C 
for  the  second  time  and  already  be  winding  directly  on  top  of  itself.  Where  and 
when  did  our  path  close?  It  closed  at  point  A,  Nqs  circuits,  before  it  reached 
point  C.  But  it  reached  point  C  (for  the  second  time)  after  NCR*NRL  circuits. 

We  see  then  that  the  exact,  integral  number  of  circuits  per  (closed)  layer  is 
given  by 

nCL  *  nCR*nRL  "  NCS 
Npc 

■  ncr<nrl  -  fi--> 


nCR‘nCS 

■  ncr(nrl-i  ♦  } 


The  exact,  nonintegral  number  of  repetitions  per  closed  layer  is  therefore  given 
by 


nRL 


Nrl-1 


*CR  ~  NCS 
*CR 
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where,  of  course,  Nrl  is  the  integral  number  of  repetitions  needed  to  close  the 
layer  and  add  N^s  additional  circuits.  Therefore,  given  N^r,  Ngo  and  Nr|_,  we 
may  compute  Ncs,  N^l ,  and  BC‘ 

We  now  show  how  to  compute  good  values  for  Nqr  and  Ngc-  Denoting  the  nomi¬ 
nal  return  angle  by  R  and  the  separation  angle  by  S 

20(L)  «  2tt*NrC  +  R 


Thus 


and 


Hence,  R  is  given  by 


Sitl 


3  NRC  +  2* 


FI -  M 

FL(--— )  =  Nrc 


R  *  2{0(L)-jt*Nrc)  =  2(0(L)-rrFL(e(L)/jr)) 

If  we  want  a  )  1  (to  avoid  the  polar  radius  lower  bound  function),  we  will  need 

R  <  S*NSC 

where 

S  *  2n/NCR 

Hence 


NSC  *  R/s  s  Ncr»R/2jt 

and  we  define  by 

Ngg  *  CE(Ncr»R/27T) 

(provided  GCF  (Ncr.Nsc)  3  1).  We  also  want  S*Nsc  -  R  to  be  as  small  as 
possible  (in  order  that  only  minimal  perturbation  of  our  nominal  quasi -geodesic 
will  be  required).  We  therefore  use  the  following  procedure  to  select  nearly 
optimal  values  of  Nqr  and  Ngc*  For  values  of  Nqr  between  some  minimum  and  maxi¬ 
mum  value,  compute: 


SO 


NSC  -  CE{NCR*R/2ff) 


'  GCF (NgQ , Nqr )  a  1,  compute 

6  *  Nsc*2jt/Ncr  -  R 

If  this  is  the  smallest  value  of  6  yet  encountered,  save  6,  NCr,  and  Ngc- 
Continue  on  to  the  next  value  of  Nqr. 

The  values  of  NCr  and  Ngc  obtained  by  this  procedure  will  insure  that  the 
return  angle  for  a  circuit  of  a  closed  repetition  will  be  as  close  as  possible 
tu  the  nominal  return  angle  while  still  maintaining  the  condition  a  >  1. 

We  now  consider  the  specification  of  a.  Our  effective  wrap  angle  function 
(t)  derivative  is  given  by 

t1 (x)  ■  aQ' (x) 

r(x)  *  /  r'(t)dt  *  /  a6'(t)dt 
o  o 


If  a  is  a  constant,  we  have 


and 


but  we  want 


t{x)  *  O0(X) 


T(L)  a  <xe(L) 


Ncr 

2T(L)  «  2n{NRC  +  sg)  +  Be 


SC,  .  2tt 


2jt(NrC  ♦  a--)  ♦  a-=-u 
RC  nCR  nCR*nCL 


Therefore 


and 


2a$iL)  •  2MBr‘  *  S'  *  *5?as 


*  *  <NRC  *  "CR  *  "ct^cl1 
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Up  to  this  point,  one  may  have  been  thinking  of  a  as  a  constant,  but  it  turns 
out  that  if  a  is  a  constant  and  r'  is  large  somewhere,  then  p'Q  will  tend  to 
follow  r'  and  be  large  also.  The  path  will  thus  deviate  considerably  from  the 
nominal  quasi -geodesic  path.  We  may  find  the  extent  of  this  deviation  intol¬ 
erable,  so  we  define  a  as  a  function: 

o(x )  *  1  ♦  iSp(x) 

where  &  is  a  constant.  We  call  p(x)  the  perturbation  function  and  we  make  p ( x > 
zero  when  we  want  pQ  to  be  exactly  equal  to  r0.  We  can  still  make  p  constant  if 
we  wish,  but  we  now  have  more  flexibility  at  a  small  computational  price.  For 
the  new  definition  of  a. 


t'(x)  «  a(x)0'(x)  *  (1+0 p(x))0'(x) 

T(x)  -  /X  T’(t)dt  »  /X  ( l-*-/Sp( t ) )0 1  ( t )dt 
o  o 

»  0(x)  +  0/X  p ( t ) 0 * (t)dt 
o 


Letting 


P(x)  »  IX  p(t)0* (t)dt 
o 

(which  can  be  computed  with  our  square  root  singularity  quadrature  formula),  we 
have 


t(x)  ■  0(x)  ♦  0P(x) 


and 


t(L)  *  0(L)  ♦  0P(l) 


so  we  can  calculate  fi  easily  as  soon  as  we  know  t(L).  Also, 

Ti  *  Ti.j  ♦  0i  -  0i-i  ♦ 

defines  the  path  data  for  our  perfectly  closed  path  in  terms 
data,  the  perturbation  function  p  and  0. 


of  our  nominal 


>ath 
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Even  though  a  is  now  a  function,  we  still  have  the  relation: 


T(LI  ■  *  safe1 


NSC  1 

ff(NRC  +  Ncr  +  NcrTNcr.Nrl-NcsI* 

Since  NRC,  Nqr,  Ngc.  and  Nc$  are  computed  on  the  basis  of  the  nominal  path,  we 
need  only  specify  NRl_  in  order  to  compute  t(L)  and  fi.  We  establish  NrL  from  the 
requirement  that  a  layer  of  filament  must  completely  cover  the  surface  for  any 
x.  We  assume  that  our  band  of  filaments  has  bandwidth  b. 


Figure  13.  Effect  of  finite  bandwidth. 
From  the  figure,  we  have  the  approximation 

b  »  c  sin  h 
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(exact  for  r  and  r0  constant;  not  valid  at  all  at  turning  points)  and 

b*  ■  c*  sin*  h  s  c*(l-cos*  h) 
but 

*o(x) 
cos  h  =  ?T5r 

Hence 

Pn  - 

b*  «  c* (1  -  (--)  ) 
but 

orr0 

Po  *  (r*;ia*-I)r*)*" 

Hence 

a  .  r*-r* 

1  -  (?9)* . 9 - 

V  '  r*  ♦  (a*-l )r* 

We  therefore  have 

c* (r*“r* ) 

hi  s 

r*+(a*-i )r* 
o 

Now,  in  order  for  one  layer  to  completely  cover  the  surface  at  this  point,  we 
clearly  must  have 

c*NCL  >  2nr 

from  which  we  conclude  that  we  need 

b,  >  ,?*!:)  — j 

Vl  'r*+(a*-l)r£' 

but  a  sufficient  condition  for  the  satisfaction  of  the  previous  inequality  is 
the  following  one: 

b*  > 

"CL  r* 

•  (i--)*(r,-r*, 

NCt  0 
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Hence 


NCl 


ro) 


but  we  want  this  inequality  to  hold  for  all  x.  Therefore,  if  we  define 

H  -  Max  (r(x) *-r0(x) * ) 

0<x<L 


then  we  need 


and 


NCt  > 


2*1^ 

‘b~“ 


nCR*nRL  •  Ncs  > 


2x/m 


We  therefore  conclude  that  a  value  of  NR^  sufficient  to  guarantee  complete 
coverage  of  the  surface  everywhere  is  given  by 


nRL 


CE((-™ 


♦  NCS)/NCR) 


Fro*  this  value  of  NRL,  we  compute  r(L)  and  fi.  If  fi  is  not  particularly  small, 
however,  we  may  be  able  to  get  away  with  an  even  smaller  value  of  NRL.  If  a 
smaller  value  of  NRL  gives  complete  coverage,  we  must  have 

a  <  0 


where 


a 


Max 

0<x<L 


P._P. 


<??-M 

NCI/  'r**(a*-l) 


-- — ) 
-l)rr 


-  b* 


The  global  maxima  M  and  m  are  easily  computed  with  a  couple  of  applications 
(global  and  local)  of  adaptive  (error  equidistributing)  sampling. 
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SURFACE  BUILDUP  RELATIONS 

In  the  previous  section,  we  discussed  how  finite  bandwidth  affects  the  sur¬ 
face  coverage  problem.  In  this  section,  we  will  discuss  how  finite  bind 
thickness  and  finite  bandwidth  affect  the  surface  buildup  problem. 

After  a  layer  of  filament  has  been  wound  onto  the  surface,  one  will  observe 
that  the  layer  varies  in  thickness  from  point  to  point.  This  is  due  to  the  fact 
that  the  hoop  angle  cannot  remain  constant  everywhere.  One  notices  in  par¬ 
ticular  that  there  is  considerable  buildup  at  and  near  the  turning  points.  The 
objective  of  this  section  is  to  compute  a  uniform  approximation  to  the  layer 
thickness  at  al_l  points  of  interest.  The  central  problem  of  this  section  is 
determining  the  length  of  a  cut,  transverse  to  the  axis  of  the  surface,  in  the 
band  of  filaments.  For  a  cylinder  and  at  a  point  where  the  hoop  angle  is 
constant,  the  length  of  such  a  cut  is  given  exactly  by 


where  b  is  the  bandwidth  and  h  is  the  hoop  angle.  At  the  turning  points,  h  is 
zero  and  the  previous  formula  predicts  infinite  buildup.  The  buildup  at  the 
turning  points  may  be  considerable,  but  it  is  certainly  not  infinite.  We  will 
leave  the  determination  of  c  for  later  and  assume  for  the  present  that  we  will 
be  able  to  obtain  it. 

We  will  first  determine  the  physical  surface  length  and  the  physical  hub 
diameters  for  a  given  bandwidth  b. 
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Figure  14.  Left  turning  point. 
Tan  a  *  r'(o)  *  ~ 

Or*  ♦  DL*  »  (|)* 
(r'(o)OL)*  ♦  DL*  a  (| ) * 


Or  » 


br‘(o)__ 

T*  r ' ( o ) * 


The  physical  length  of  the  surface  is  therefore  given  by 


l  ♦  \  (-— i--r-  ♦ 

/l+r'(o)*  /l+r’ (L)a 


and  the  left  and  right  physical  hub  diameters  are  given  by 


and 


2r(o) 


_br^o] _ 

/l+r* (o)* 


2r(L)  ♦ 


..brlltl. 

/l*r'(L)« 


respectively. 

In  what  follows,  we  will  let  the  6  quantities  refer  to  a  single  band  of 
filaments  and  the  A  quantities  refer  to  the  aggregation  of  these  bands  which 
form  the  layer. 

From  the  same  diagram  we  also  have 

t  »  ir  cos  a 


so  that 

5r  *  t  /l+r ' (x)* 

for  any  point  x  and  any  orientation  of  the  band. 

We  make  the  following  assumption  or  idealization  in  our  modeling  of  the 
behavior  of  the  band  of  filaments.  We  assume  that  the  band  behaves  as  a 
flexible* membrane  of  constant  thickness  t  and  width  b  which  clings  uniformly  to 
the  surface.  The  area  of  a  transverse  cut  through  the  band  is  therefore  given 
by 

(ir(r+5r )  *-nr* )  *  6A 


where  6y  is  the  angular  length  of  the  cut.  Since  c  *  r6y,  we  have 

~  (2r6r+5r* )  *  «A 
zr 
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/ 


We  Hill  now  have  two  Aa  area  contributions  from  each  circuit,  hence,  the 
total  area  contribution  for  our  layers  of  wrap  is  given  by 

AA  *  2N£lAA 

Now,  if  we  think  of  Ar  as  the  average  thickness  of  aggregate  filament  material 
buildup  at  any  point,  we  have 

AA  *  »r(r+Ar)*  -  er* 

■  Jt(2rAr+Ar*) 

Solving  this  quadratic  equation  for  Ar  and  rational izing,  we  have 


Ar .  ...ays . 

r  /r*+AA/7r 

Technically,  the  previous  value  of  AA  is  due  to  filaments  alone,  while  we  should 
also  consider  the  matrix  material  component.  If  we  let  v  be  the  fi lament /volume 
ratio,  then  AA  for  matrix  and  filament  together  is  given  by 

AA  »  2NclAA/i> 


In  summary,  we  have 


A  r  a 


AA/7T 


r  ♦  /r*+AA/ff 
AA  ■  2NCL«5A/V 
AA  «  cAr(l+Ar/{2r) ) 


Ar  *  t/l+r ' 1 


The  only  quantity  in  these  equations  whose  computation  has  not  yet  been 
addressed  is  the  length  of  the  transverse  filamentary  cut,  c.  We  will  now 
proceed  to  compute  c.  In  order  to  complete  our  model  of  the  band,  we  need  to 


define  the  family  of  paths  orthogonal  to  our  quasi -geodesic.  These  orthogonal 


paths  will  enable  us  to  define  precisely  the  finite  width  of  our  band  and  com- 


We  let  s  denote  the  length  of  our  quasi -geodesic  path,  a  be  the  length  of  the 
corresponding  orthogonal  path,  and  m  be  the  length  of  the  meridian.  The  hoop 
angle  H  of  the  orthogonal  path  is  defined  by 

cos  h  =  sin  H 


6i 


V 


The  contribution  to  the  meridian  is 


but 


dm  a  da  sin  H  =  do  cos  h  * 


r 


da 


Therefore 


dm*  »  dx*  >  dr*  a  (l+r'*)dx* 


da 


dm  *  —  M+r '  *  dx 


hence 

c»*  «  --  /i+rv"* 

ro 

We  call  the  orthogonal  path  wrap  angle  n  and  we  have 

-rdrj  a  da  cos  H  «  da  sin  h 
*  do(l-cos*  h}3* 

«  da(l-(r0/r)*)^ 

*  f  da(r*-r*),< 

i  r  o 

Hence  | 

dr>  =  -  J;  (r^r*)3*  •  £-  (l+r1  *)%  dx 

and 

!  H'  *  -  {(r*-r*)(l4r'*)}J| 

Note  that  n  is  a  decreasing  function  of  x. 


6! 


point  on  axis 


orthogonal 

paths 


Figure  16.  Determination  of  cut  length. 

Going  from  A  to  E  indirectly  by  the  route  ABCOE,  we  pass  through  the  wrap  angle 
6y.  which  is  decomposed  in  terms  of  a  and  s  wrap  angles: 

6y  •  6^1  +  +  602  ♦  0>12 

X-0X1  x  X+0X2  x 

*  /  17'  (t)dt  +  /  0'  (t)dt  ♦  /  0*(t)dt  +  /  o’(t)dt 

x  X-6X1  x  X+0X2 

x+0X2 

-  /  0'(t)  -  r?'(t)dt 

x-6xj 

where  x  is  the  axial  location  of  the  cut  and  6x1  and  6x2  are  the  axial  distances 
therefrom  to  points  B  and  D,  respectively.  Now 


••  -*■  -  ? 

0 

♦  { (r*-r*)  (1+r ' * ) }% 


1  r«  ( r  *  —r  *  J5* 

l  (1+P..,*{  2  + 

r  ( r- 2 -Pq  )  ^  ro 


r*  ♦  r*  -  p* 
o  o 


-  I  (l+r’ * }*< 

r0(r*-r*)^ 


C.  ,lt£ :l!)H  ,  r. 

o  1  1 Q 


(note  that  the  only  difference  between  y’  and  8’  is  that  0’  contains  the  factor 
r0/r  and  y 1  contains  the  factor  r/r0).  The  length  of  the  transverse  cut  in  the 
band  is  therefore  given  by 

x+6x  2 

c(x)  *  r(x)6y  *  r(x)  /  y’(t)dt 

x-Sxj 


where 


*  r(x)(y(x+<5x2)-y(x-6x1)) 

y* (x,  *  cl? i. 

T  ' '  rft(x)  'r(x)l-rn(x)* 


and  5xj  and  5x2  are  determined  by 


,x  b  ,*+«*2 

/  o' (t)dt  =  f  *  / 
x-6xi  2  x 


<7’(t)dt 


where 


(l*r'(x)*)* 


In  the  special  case  when  r  and  r0  are  constant,  we  have  as  a  check 


-  r  b  ,  r 
1  ro  2  z  r0 

bro 

5xj  x  6x2  = 

y- . C . . 1 . 

ro(r*-rdH  r0(l-cos*h)* 


rQ  sin  h 


=  const 
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i 


X+dX£ 

c(x)  «  r  /  y'(t)dt  *  r*26xi«y' 
x-6xj 


»  2r 


_ 1 _ 

r0  sin  h 


__b__ 
*  sin  h 


Since  o'  is  always  greater  than  unity,  a  is  strictly  monotone  increasing  end 
hence  has  an  inverse  a-1.  We  can  therefore  compute  fixj  and  5x2  in  the  following 
manner: 

<r(x)  -  aix-Sx^  «  | 

<y(x-5xj)  *  o(x)  -  | 
x  -  5xj  •  -  |) 

5xx  «  x  -  o~'(o(x)  -  |) 

Also 

o(x+ 6x2)  -  ff(x)  »  | 

o(x+6x2)  =  <7{X)  ♦  | 

x  ♦  «x2  *  o~'(o(x)  *  |) 

i  «X2  *  <J-'[o[x)  +  |)  -  X 

If  we  approximate  a’  by  a  piecewise  linear  function,  we  have  another  application 
for  our  process  of  inverting  the  integral  of  a  positive,  continuous,  piecewise 
linear  function.  We  only  need  to  decide  what  mesh  to  use  for  a'. 

Let! hatted  items  denote  piecewise  linear  estimates  and  unhatted  items 
denote  e^act  values.  The  error  in  computing  6xj  is  therefore: 
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Letting 


e(x)  »  5xj  -  5xj 
«  x  -  o~ 1  (o(x )  -  \) 

-  (x-ff-’(o{x)  -  |)) 

■  5"’(o(x)  -  -  a_,(ff(x)  -  |) 


we  have 


o' 1  ■  S  and  a-1  »  S 


e(x)  «  S(a(x)  -  -  S (a(x)  - 

Assuming  b  to  be  small,  we  have 

*(x)  «  S(o(x) )  -  |  S' (o(x) ) 
-  S(o(x) )  ♦  |  S’  (<r(x) ) 
but 


S(0(X) )  *  X  *  S(a(x) ) 

hence 

*(x)  »  |  (S' (o(x) )-§’ (o(x) ) ) 

Now,  since 

ff'Mx)  »  S(X)  , 

X  a  <JT(S(X)  ) 

therefore 

1  «  (o*(S(x))S'(x) 
or 

I 

S' (x)  ±  l/o' (S(x) ) 

1 

Also 

\ 

S'(0(X))  a  I/O' (S(o(x) ) )  a  l/0'(x) 
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rot  O' 


Similarly 


A  rough  bound  for  e  on  the  i*^  sub-interval  is  therefore 

I  •  I  <  |c'i  •  l  hil»’”l(i) 

In  order  to  make  e  roughly  constant,  we  choose  our  g  function  as 

g(x)  «  j  a"  1  (x)  |  Vo'  (x) 

Let  us  now  obtain  an  asymptotic  expression  for  the  length  of  the  cut  at  the 
left  turning  point  (for  one-half  circuit).  From  geometric  considerations  we 
know  that  5xj  *  0.  We  approximate  6x2  using  a  two-term  Taylor  series: 

5x2  -  S(o(x)  +  |)  -  x 

-  S(o(x) )  ♦  S'  (cr(x)  >  |  -  x 
2o”(x) 


We  use  the  previously  derived  quadrature  rule  to  integrate  y': 


c(o)  »  r(o)  /  y’(x)dx  «  r(o)  /  2 

o  o  ^ 

a  r(o)  •  2/6x2  /6x2/3  y'(5x2/.3) 


2r(o)5x2 

~/F“" 


y' (6x2/3) 
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but 
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If  r'(o)  is  large,  we  have 


c(o)  »  (br(o))H 

If  we  use  virtually  the  same  analysis  to  approximate  the  cut  at  a  distance  of 
b/2  from  the  left  turning  point  meridionally,  we  get 

c(x)  -  fi  c(o) 

The  length  of  the  cut  v*.  therefore  about  40  percent  greater  here  than  at  the 
turning  point. 

FUNCTION  DEFINITION 

The  off-line  filament  winding  process  begins  with  the  explicit  specifica¬ 
tion  of  the  profile  or  radius  function  of  the  surface  of  revolution  to  be  wound. 
We  have  been  calling  this  function  r(x).  Although  one  can,  by  exercising  proper 
caution,  wind  across  discontinuities  in  r'  and  r,  we  require  here  for  the  sake 
of  simplicity  that  r '  exist  everywhere  in  (0,L).  This  is  not  to  say,  however, 
that  we  cannot  allow  r*  and  r"  to  be  fairly  large  in  some  spots;  we  can. 

The  next  step  is  to  define  the  desired  nominal  path  we  wish  the  filament  to 
follow  on  the  surface.  This  is  done  by  explicitly  specifying  the  polar  radius 
function  r0(x).  For  any  acceptable  pure  geodesic  path,  r0(x)  is  of  course  just 
a  constant.  For  a  quasi -geodesic  path,  however,  we  may  need  as  much  flexibility 
defining  r0  as  we  need  defining  r. 

The  last  function  we  need  to  specify,  the  perturbation  function  p(x),  is 
required  for  flexibility  in  specifying  the  slight  modification  which  we  must 
make  in  our  nominal  quasi -geodesic  path  in  order  to  wrap  the  surface  with  a 
closed,  uniformly  covering  layer  of  filament.  In  many  cases,  this  function 
could  simply  be  defined  as  a  constant,  but  in  the  event  that  r'  is  large  near 
the  turning  points,  a  constant  pwill  have  the  undesirable  effect  of  forcing  the 
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effective  polar  radius  function  pQ(x)  to  mimic  the  behavior  of  r  near  the 
turning  points.  Making  p  zero  in  a  neighborhood  of  the  turning  points  leaves 
pQ  equal  to  r0  in  these  neighborhoods. 

The  smoothness  of  a  function  is  specified  in  the  following  Banner.  A  func¬ 
tion  f  is  said  to  be  Cn  if  at  least  the  first  n  derivatives  of  f  exist  and  are 
continuous.  The  most  stringent  requirement  on  the  smoothness  needed  by  r,  r0, 
and  p  can  be  determined  by  examining  the  g  function  for  generating  the  winder 
data.  The  extent  of  differentiation  in  this  function  indicates  that  r"',  rg, 
and  p”  should  at  least  exist  almost  everywhere;  r  should  be  C*  and  r0  and  p 
should  be  C1.  Additionally,  if  we  want  g  to  be  continuous,  we  will  need  r  to  be 
C*  and  r0  and  p  to  be  C*. 

Of  course  it  is  not  at  all  necessary  for  actual  machining  of  the  surface  to 
be  done  this  smoothly,  but  it  is  necessary  for  the  functions  we  compute  with  to 
be  this  smooth  in  order  to  give  us  computationally  reliable  algorithms. 

There  are  any  number  of  ways  in  which  we  might  define  these  functions.  The 
method  chosen  here  is  to  begin  with  a  piecewise  linear  function  and  round  off 
the  corners  to  obtain  the  desired  smoothness.  We  do  the  rounding  off  by  using 
smoothing  by  averaging  (refs  3,4).  Consider  the  smoothing  operator  S  defined  by 

1  »x+h 

Sf(x)  •  \ -  f  f(t)dt  h  >  0 
2h  x-h 

Differentiating  this  equation  with  respect  to  x  using  Leibnitz’s  rule  givas  us: 

3x  Sf(x)  *  !h  <f<*+h)-f(x-h)) 

3H.  S.  Shapiro,  Smoothing  and  Approximation  of  Functions,  Van  Nostrand  Reinhold 
Company,  New  York,  1969. 

4R.  W.  Soanes,  "Function  Smoothing  by  Repeated  Averaging,"  ARDEC  Technical 
Report  ARCCB-TR-88012,  Benet  Laboratories,  Watervliet,  NY,  March  1988. 
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and 


_d' 

dx1 


n+1 

— -  Sf(x) 
n+1 


Hence,  we  see  that  if  f  is  Cn,  then  Sf  is  Cn+1.  If  we  apply  S  i  times,  we  also 
have  that  if  f  is  Cn,  then  S’f  is 

We  use  the  following  symbolism  for  the  successive  integrals  of  a  function 
f: 


f0<*)  ■  f(x) 

f i (x)  -  r  f i-! (t )dt  i  >  1 

c 

By  definition  of  S, 

Sf (x)  *  /***  f(t)dt  »  <fl(x*h)-fl(x-h» 

1  ,**h 

S*f (x)  «  Jr  /  Sf(t)dt 
2h  x-h 

» k  r_l  1r  (fift^j-f^t-hiidt 

i  ,x+h  ,x*h 

■  ih  lfx.h  'l'*-""") 

i  ,x+2h  ,x 

-  iRi  (fx  fH*>«  -  /x.2hfl(t)dt) 

*  iRI  (f2(x*2h)-f2(x)-(f2(x)-f2(x-h))) 

»  4RI  {f2(x+2h)-2f2(x)+f2(x-2h) ) 


In  general,  it  can  be  proved  (ref  4)  that 

i 

Sif(x)  =  -™-r  (-l)k(,J)f  i  (x+(i-2k)h) 

k=0 


*£7  W.  Soanes,  "Function  Smoothing  by  Repeated  Averaging,"  ARDEC  Technical 
Report  ARCCB-TR-88012,  Benet  Laboratories,  Watervliet,  NY,  March  1988. 
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For  the  rounded  area  in  the  vicinity  of  X(<,  we  will  take  c  *  xj<.  This 
establishes  a,b  and  f-j  for  any  value  of  i.  In  order  to  subsequently  specify 
S’f,  we  need  to  establish  the  one  remaining  parameter  h.  A  point  x  will  be 
called  a  linear  point  of  S’f  if  x+ih  and  x-ih  both  lio  in  the  domain  of  the  same 
linear  section  of  the  initial  piecewise  linear  approximation.  We  refer  to  such 
a  point  as  linear  because  S’  preserves  linear  functions  exactly,  but  since  our 
basic  approximation  is  only  piecewise  linear.  S’  will  only  preserve  the  function 
at  the  point  x  if  x*ih  and  x-ih  both  lie  in  the  domain  of  the  same  linear  sec¬ 
tion. 

We  want  at  least  the  midpoints  of  each  subinterval  to  be  linear  points  so 
that  we  may  have  different  values  of  h  associated  with  each  node  or  rounded  sec¬ 
tion.  We  therefore  want  the  following  two  inequalities  to  be  satisfied: 
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This  leads  us  to  define  the  h  value  associated  with  xk  to  be 

h  -  ft  K  -  xk.t  ,  xkn  -  xk} 

or  some  fraction  thereof.  Hence,  S’f  is  defined  between  the  Midpoints  of  the 
subintervals  surrounding  xk. 

INDEFINITE  INTEGRATION 

In  this  section  we  define  and  analyze  the  complexity  of  an  algorithm  for 
indefinite  numerical  integration  which  tends  to  minimize  the  effects  of  round¬ 
off  error  in  the  limit  as  the  number  of  mesh  points  becomes  large.  Typically, 
we  want 

F(x)  -  J*  f <t)dt 
*1 

estimated  at  the  mesh  points  xj,  xj,  ...xn,  where  n  can  be  large.  This  can  be 
accomplished  by  the  usual  recursion. 

,*i 

Fi  *  Fi-1  +  /  f(x)dx  i  «  2,3,... ,n 
xi-l 

where  Fj  «  0  and  we  estimate 

*i 

/  f (x)dx 
*i-l 

by  any  desired  quadrature  rule. 

If  n  is  large,  the  usual  recursion  will  suffer  from  round-off  error  as  we 
continually  add  relatively  small  quantities  to  increasingly  larger  ones.  The 
larger  n  is  and  the  greater  the  value  of  the  accumulated  integral,  the  greater 
will  be  the  number  of  significant  digits  effectively  dropped  from  each  subinter 
val's  contribution  to  the  integral.  It  should  be  clear  that  Fn  will  suffer  the 
most  from  round-off.  To  reduce  the  effects  of  round-off  error,  we  construct  an 
algorithm  which  will  tend  to  add  only  quantities  of  similar  magnitude  and 
thereby  make  use  of  roughly  the  same  number  of  significant  digits  in  each 
addend.  Consider  the  following  algorithm: 
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T 


Initiation  step: 

F(  1 ) : *0 

i :  *2 

do  until  i  >  n 
{L{ i > :«i-l 
*i 

F ( i ):*/  f(x)dx 

xi-l 

i:«i*l} 

Summation  step: 

do  until  l(n)*l 
{ i  :*n 

do  until  L(i)=l 
(k:*L( i ) 

F( i )  :*F( i )  ♦  F\k) 

L ( i ) :*L(k) 
i :«i-l } } 

Examining  the  summation  step,  we  see  that  during  each  pass  of  this  step  we 
have  that  F(i)  and  F(k)  each  represent  integrals  over  the  same  number  of  subin¬ 
tervals,  providing  L(k)  *  1.  Hence,  the  new  value  of  F(i)  will  be  the  sum  of 
two  numbers  of  similar  magnitude  (except  for  the  last  time  F(i)  is  updated).  We 
now  examine  the  space-time  complexity  of  the  usual  recursion  and  the  new 
algorithm.  The  usual  recursion  needs  an  array  of  size  n  and  the  new  algorithm 
needs  two  arrays  of  size  n.  Hence,  both  these  algorithms  have  the  same  space 
complexity,  0(n).  The  number  of  floating-point  additions  needed  in  the  usual 
recursion  is  n-2;  hence,  the  time  complexity  of  the  usual  recursion  is  0(n).  We 
will  now  analyze  the  time  complexity  of  the  new  algorithm. 

We  first  note  how  L  changes  from  pass  to  pass  via  Table  I  (n=20). 
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TABLE  1.  BEHAVIOR  OF  L  ARRAY  FROM  PASS  TO  PASS 


x  indices 

L  initially  : 
L  after  first  pass  : 
L  after  second  pass: 
L  after  third  pass  : 
L  after  fourth  pass: 
L  after  fifth  pass  : 


On  the  basis  of  this  table,  we  can  construct  the  following  table. 

TABLE  II.  AOOITIONS  PER  PASS 
p  (pass  index)  ap  (floating-point  adds  per  pass) 


1 

2 

3 

4 

5 

• 

P 


n-2 

n-2-1 

n-2-1-2 

n-2-1-2-4 

n-2-1-2-4-8 

• 

n-2-1-2-4-... -2P‘2 


Me  therefore  have  that 


by  induction,  but 


P-2 

ap  *  "-2-  J  2’ 

i*0 

k 

S*  «  J  2’  ■  2k+l  -  1 

i*0 


% 
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Therefore 


ap  *  n-2-Sp_2 

*  n-2- ( 2p— 1 - 1 ) 


In  all  that  follows,  log  will  denote  logarithms  to  the  base  2  ard  In  will  denote 
logarithms  to  the  base  e.  For  each  pass,  including  the  last,  we  have 

•p  >  0 

n-l-2p_1  >  0 
2p-1  <  n-1 
p-1  <  log(n-l) 

j' 

p  <  1  ♦  log(n-l) 

The  cumulative  number  of  adds  in  p  passes  is 

P  P 

Ap  -  ^  ai  *  J  (n-1-21'1) 
i*l  i*l 


P 

«  p(n-l)  -  J  2i_1 
i*l 


p-1 

*  P(n-l)  -  ])  21  «  p(n-l)  -  Sp.j 

i*0 


*  p(n-l)  -  (2P-1) 

The  total  number  of  f 1oating-point  additions  in  the  summation  step  is  therefore 

Ap  *  p(n-ll-2P+l 

where  p  is  the  largest  inteqer  strictly  less  than 

1  ♦  log(n-l) 
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If  n-1  is  an  exact  power  of  2, 


p  *  log(n-l) 
and 

Ap  »  (n-l)log(n- 1)  -  (n-1)  +  1 
■  (n-l)(log(n-l)-log  2)  +  1 
«  (n-l)log(-”)  ♦  1 
If  n-1  is  not  an  exact  power  of  2, 

p  *  1  ♦  log(n-l)  -  e 

where  e  is  some  positive  fraction.  Hence 

Ap  *  (n-1) (l+log(n-l)-e) 

_  2l*1og(n-l)-e  +  1 

•  (n-1) (l+1og(n-l)-e) 
-21_e(n-l )  ♦  1 

-  (n-1) (l+log(n-l)-e-21~€)  ♦  1 
To  get  an  upper  bound  on  Ap,  we  need  an  upper  bound  on 

f ( £ )  =  -c-21"6 

Note  that  f ( 0 )  *  f ( 1 )  »  -  2.  Differentiating,  we  have 

f’(£)  *  -1  -  (-l)2W1n2 


«  -1  + 


lilt. 

Tog  e 


and 


f"(€) 


log  e  (Tog  e) ‘ 


Since  f"(£)  <  0  for  all  e,  f(e)  attains  a  maximum.  Setting  f'(e)  »  0  gives  us 


21_€  *  log  e 
1-e  *  log  log  e 
c  *  1  -  log  log  e 
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v--  1 


hence 


f(e)  =  -1  +  log  log  e  -  log  e 
*  -1  +  log 

Since 

Ap  ■  (n-l)(l+1og(n-l)+f (€))  ♦  1 
for  the  correct  e,  we  have 

Ap  <  (r.-l)  (l+1og(n-l)-l+ log  )  +  1 

for  the  maximizing  e.  Therefore, 

Ap  <  {n-l)log(j^)  ♦  1 

<  (n-l)log(j"i|j)  .  1 

which  is  our  sought  upper  bound  on  the  time  complexity  of  the  summation  step. 

We  see  that  the  price  we  have  to  pay  for  the  reduction  of  round-off  error 
is  the  logarithmic  factor  in  our  bound.  Note,  however,  that  this  factor  grows 
very  slowly  with  n. 
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220  7TH  STREET,  N.E. 

CHARLOTTESVILLE,  VA  22901 

COMMANDER 

US  ARMY  LABCOM 

MATERIALS  TECHNOLOGY  LAB 

ATTN:  SLCMT-IML  (TECH  LIB)  2 

WATERTOWN,  MA  02172-0001 


* 


NOTE:  PLEASE  NOTIFY  COMMANDER,  ARMAMENT  RESEARCH,  DEVELOPMENT,  AND  ENGINEERING 
CENTER,  US  ARMY  AMCCOM,  ATTN:  8ENET  LABORATORIES,  SMCAR-CCB-TL , 
WATERVLIET,  NY  12189-4050,  OF  ANY  ADDRESS  CHANGES. 


TECHNICAL  REPORT  EXTERNAL  DISTRIBUTION  LIST  (CONT'O) 


NO.  OF 
COPIES 

COMMANDER 

US  ARMY  LABCOM,  ISA 

ATTN:  SLCIS-IM-TL  1 

2800  POWDER  MILL  ROAD 
ADELPHI,  MD  20783-1145 

COMMANDER 

US  ARMY  RESEARCH  OFFICE 

ATTN:  CHIEF,  IPO  1 

P.O.  BOX  12211 

RESEARCH  TRIANGLE  PARK,  NC  27709-2211 
DIRECTOR 

US  NAVAL  RESEARCH  LAB 
ATTN:  MATERIALS  SCI  &  TECH  DIVISION 
CODE  26-27  (DOC  LIB) 

WASHINGTON,  O.C.  20375 


NO.  OF 
COPIES 

COMMANDER 

AIR  FORCE  ARMAMENT  LABORATORY 
ATTN:  AFATL/KN  1 

EGLIN  AFB,  FL  32542-5434 

COMMANDER 

AIR  FORCE  ARMAMENT  LABORATORY 
ATTN:  AFATL/MNF 

EGLIN  AFB,  FL  32542-5434  1 

METALS  AND  CERAMICS  INFO  CTR 
BATTELLE  COLUMBUS  DIVISION 
505  KING  AVENUE 

COLUMBUS,  OH  43201-2693  1 


NOTE:  PLEASE  NOTIFY  COMMANDER,  ARMAMENT  RESEARCH,  DEVELOPMENT,  AND  ENGINEERING 
CENTER,  US  ARMY  AMCCOM,  ATTN:  BENET  LABORATORIES,  SMCAR-CCB-TL, 
WATERVLIET,  NY  12189-4050,  OF  ANY  ADDRESS  CHANGES. 


